On the Numerical Solution of the Exact Factorization Equations

TL;DR

This paper demonstrates a self-consistent numerical solution to the exact factorization equations, revealing their stability properties and implications for non-adiabatic molecular dynamics simulations.

Contribution

It introduces a novel approach to solve the EF equations self-consistently, addressing stability issues due to their non-Hermitian form.

Findings

Stable propagation achieved in model systems

Identification of non-Hermitian Hamiltonian effects

Insights into non-adiabatic dynamics modeling

Abstract

The exact factorization (EF) approach to coupled electron-ion dynamics recasts the time-dependent molecular Schr\"odinger equation as two coupled equations, one for the nuclear wavefunction and one for the conditional electronic wavefunction. The potentials appearing in these equations have provided insight into non-adiabatic processes, and new practical non-adiabatic dynamics methods have been formulated starting from these equations. Here we provide a first demonstration of a self-consistent solution of the exact equations, with a preliminary analysis of their stability and convergence properties. The equations have an unprecedented mathematical form, involving a non-Hermitian Hamiltonian, and so the usual numerical methods for time-dependent Schr\"odinger fail when applied in a straightforward way to the EF equations. We find an approach that enables stable propagation long enough to witness non-adiabatic behavior in a model system before non-trivial instabilities take over. Implications for the development and analysis of EF-based methods are discussed.

Figures

Equations

$$\Psi (r,R,t)={\varphi }_R(r,t)\chi (R,t)$$

$$\int {\left|{\varphi }_R(r,t)\right|}^{2}dr=1,\forall R,t.$$

$${\varphi }_R(r,t)$$

$$\chi (R,t)$$

$$\left[{\hat{H}}_{\mathrm{BO}}+{\hat{V}}_{\mathrm{ext}}^e(r,t)+{\hat{U}}_{\mathrm{en}}-ϵ(R,t)\right]{\varphi }_R(r,t)=i{\partial }_t{\varphi }_R(r,t),$$

$$\left[\frac{{\left(-i{\nabla }_R+A(R,t)\right)}^2}{2M}+V_{\mathrm{ext}}^n(R,t)+ϵ(R,t)\right]\chi (R,t)=i{\partial }_t\chi (R,t),$$

$$ϵ(R,t)$$

$$A(R,t)$$

$${\hat{U}}_{\mathrm{en}}[\varphi ,\chi ]$$

$$A(R,t)\to A(R,t)+{\nabla }_R\theta (R,t)$$

$$ϵ(R,t)\to ϵ(R,t)+{\partial }_t\theta (R,t)$$

$$0\stackrel{\mathrm{PNC}}{=}\frac{\partial }{\partial \lambda }⟨{\varphi }_R(t)|{\varphi }_R(t){⟩}_r=-\mathrm{Im}{⟨{\varphi }_R(t)|-i{\partial }_{\lambda }{\varphi }_R(t)⟩}_r$$

$$⟨{\varphi }_{R+\Delta R}(t)|{\varphi }_{R+\Delta R}(t)⟩-⟨{\varphi }_R(t)|{\varphi }_R(t)⟩\to$$

$$\frac{1}{\Delta R}\sum _i|{\varphi }_{i,J+1}^{(n)}{|}^2-|{\varphi }_{i,J}^{(n)}{|}^2$$

$$=2\mathrm{Re}\left(\sum _i\frac{({\varphi }_{i,J+1}^{(n)\ast }+{\varphi }_{i,J}^{(n)\ast })}{2}\frac{({\varphi }_{i,J+1}^{(n)}-{\varphi }_{i,J}^{(n)})}{\Delta R}\right)$$

$$=-2\mathrm{Im}\left(\sum _i\frac{({\varphi }_{i,J+1}^{(n)\ast }+{\varphi }_{i,J}^{(n)\ast })}{2}\frac{-i({\varphi }_{i,J+1}^{(n)}-{\varphi }_{i,J}^{(n)})}{\Delta R}\right)$$

$$\begin{gathered} ⟨{\varphi }_R(t+\Delta t)|{\varphi }_R(t+\Delta t)⟩-⟨{\varphi }_R(t)|{\varphi }_R(t)⟩\to \\ \sum _i\left(|{\varphi }_{i,J}^{(n+1)}{|}^2-|{\varphi }_{i,J}^{(n)}{|}^2\right)/\Delta t \\ =2\mathrm{Im}\left(-i\sum _i\frac{({\varphi }_{i,J}^{(n+1)\ast }+{\varphi }_{i,J}^{(n)\ast })}{2}\frac{({\varphi }_{i,J}^{(n+1)}-{\varphi }_{i,J}^{(n)})}{\Delta t}\right) \end{gathered}$$

$$\begin{gathered} \sum _i\frac{({\psi }_i^{(n+1)\ast }+{\psi }_i^{(n)\ast })}{2}\frac{-i({\psi }_i^{(n+1)}-{\psi }_i^{(n)})}{\Delta t} \\ =\sum _{i,l}\frac{({\psi }_i^{(n+1)\ast }+{\psi }_i^{(n)\ast })}{2}{H}_{i,l}\frac{({\psi }_l^{(n+1)}+{\psi }_l^{(n)})}{2} \end{gathered}$$

$$\begin{gathered} \left({\hat{H}}_{\mathrm{BO}}+V_{\mathrm{ext}}^e\right){\varphi }_R(r,t)-⟨{\varphi }_R(t)|{\hat{H}}_{\mathrm{BO}}+V_{\mathrm{ext}}^e|{\varphi }_R(t){⟩}_r{\varphi }_R(r,t) \\ -\frac{1}{2M}{\nabla }_R^2{\varphi }_R(r,t)+⟨{\varphi }_R(t)|\frac{1}{2M}{\nabla }_R^2|{\varphi }_R(t){⟩}_r{\varphi }_R(r,t) \\ -\frac{1}{M}\frac{{\nabla }_R\chi (R,t)}{\chi (R,t)}\cdot {\nabla }_R{\varphi }_R(r,t)+i\left[\frac{1}{M}\frac{{\nabla }_R\chi (R,t)}{\chi (R,t)}\cdot A(R,t)\right]{\varphi }_R(r,t) \\ =i{\partial }_t{\varphi }_R(r,t)-⟨{\varphi }_R(t)|i{\partial }_t|{\varphi }_R(t){⟩}_r{\varphi }_R(r,t), \end{gathered}$$

$$|{\varphi }_R(t)⟩\equiv \frac{({\varphi }_{i,J}^{(n+1)}+{\varphi }_{i,J}^{(n)})}{2}$$

$$\begin{gathered} \tilde{A}(R,t)\equiv \frac{-i}{{\mathcal{N}}^2(R,t)}\times \\ \sum _i\frac{({\varphi }_{i,J}^{(n+1)}+{\varphi }_{i,J}^{(n)})}{2}\frac{({\varphi }_{i,J+1}^{(n+1)}-{\varphi }_{i,J}^{(n+1)})+({\varphi }_{i,J+1}^{(n)}-{\varphi }_{i,J}^{(n)})}{2\Delta R} \end{gathered}$$

$${\mathcal{N}}^2(R,t)=\sum _i\frac{({\varphi }_{i,J}^{(n+1)\ast }+{\varphi }_{i,J}^{(n)\ast })}{2}\frac{({\varphi }_{i,J}^{(n+1)}+{\varphi }_{i,J}^{(n)})}{2}$$

$${\hat{U}}_{\mathrm{en}}\to {\hat{U}}_{\mathrm{en}}-i\mathrm{Im}⟨{\Phi }_R(t)|{\hat{U}}_{\mathrm{en}}|{\Phi }_R(t)⟩.$$

$${\varphi }_R(r,t)=\sum _j{C}_j(R,t){\varphi }_R^{j,BO}(r)$$

$$i{\partial }_t{C}_i(R,t)$$

$$U_{\mathrm{en}}^i$$

$$d_{ij}=d_{ij}(R)=⟨{\varphi }_R^i|{\nabla }_R{\varphi }_R^j⟩,d_{ij}^2=d_{ij}^2(R)=⟨{\varphi }_R^i|{\nabla }_R^2{\varphi }_R^j⟩$$

$$A(R,t)$$

$$ϵ(R,t)$$

$$⟨{\varphi }_R|{\partial }_t{\varphi }_R⟩$$

$$⟹ϵ(R,t)$$

$${\partial }_t{C}_i(R,t)\to {\partial }_t{C}_i(R,t)-⟨{\varphi }_R|{\partial }_t{\varphi }_R{⟩}_r{C}_i(R,t)$$

$$H(r,R)$$

$$-\frac{\mathrm{erf}\left(\frac{|R-r|}{R_f}\right)}{|R-r|}-\frac{\mathrm{erf}\left(\frac{|r-\frac{L}{2}|}{R_r}\right)}{|r-\frac{L}{2}|}-\frac{\mathrm{erf}\left(\frac{|r+\frac{L}{2}|}{R_l}\right)}{|r+\frac{L}{2}|}$$

$$i{\partial }_t\vec{C}$$

Full text

I Introduction

The exact factorization (EF) of the time-dependent molecular Schrödinger equation Abedi et al. (2010, 2012) is an exact reformulation of the quantum dynamics of interacting electronic and nuclear systems. The molecular wavefunction is expressed as a single product of a nuclear wavefunction and an electronic wavefunction that is conditionally dependent on the nuclear coordinate, with two coupled equations describing their motion. Many interesting exact properties of the formalism have been uncovered in recent works, e.g. Abedi et al. (2010, 2012); Suzuki et al. (2014); Khosravi et al. (2015); Abedi et al. (2013a); Min et al. (2014); Curchod et al. (2016); Fiedlschuster et al. (2017); Requist et al. (2016, 2017); Curchod and Agostini (2017); Lefebvre (2015); Agostini and Curchod (2018); Scherrer et al. (2017); Eich and Agostini (2016); Schild et al. (2016); Requist and Gross (2016); Li et al. (2018); Cederbaum (2015); Hoffmann et al. (2018); Schild and Gross (2017); Gonze et al. (2018), shedding light on the nature of interactions between dynamical quantum subsystems as well as on interactions between quantum and classical subsystems beyond adiabatic treatments. From a practical viewpoint, the EF equations provide a rigorous starting point for methods for non-adiabatic dynamics and already we have seen the development of mixed quantum-classical approaches Agostini et al. (2016); Min et al. (2015); Gossel et al. (2018); Agostini et al. (2018), with successful applications in photochemical dynamics Min et al. (2017), as well as density functionalizations Requist and Gross (2016); Li et al. (2018).

Regarding the exact EF equations, prior studies of the exact features were based on first solving the original full molecular Schrödinger equation, and then extracting the exact coupling terms from inverting the exact EF equations. That is, a direct numerical solution of the coupled exact EF equations was avoided. Indeed, the stability and convergence properties of these equations remained unexplored, properties which are of interest when developing further EF-based approximations. In this paper, we discuss unique challenges that a self-consistent numerical solution of the exact coupled EF equations pose, give a numerical solution for a model problem, and present a preliminary analysis of their stability. We show that the usual numerical methods developed for the time-dependent Schrödinger equation fail when applied to the EF equations, and investigate these failures with a preliminary formal analysis. We show how we were able to obtain a stable numerical propagation for long enough to witness non-adiabatic behavior in the model system before non-trivial instabilities kill the calculation.

This paper is structured as follows: section II briefly reviews the EF formalism. Section III then discusses aspects that a numerical solution need to consider, before discussing such a solution for a model problem in Sec. IV. A preliminary analysis of the well-posedness of the EF equations is given in Sec. V, before we conclude in Sec. VI.

II The exact factorization equations

In the EF approach Abedi et al. (2010, 2012); Gidopoulos and Gross (2014); Hunter (1975, 1981), the exact molecular wavefunction is written as a single product of a marginal wavefunction $\chi$ and a conditional wavefunction $\Phi_{R}$, as

[TABLE]

subject to the partial normalization condition (PNC)

[TABLE]

Here, $R$ represent the set of all nuclear coordinates and $r$ the electronic coordinates. The factorization above is unique up to phase-transformation that depends only the nuclear coordinates and time:

[TABLE]

Notice that Eq. (1) has the same form as the wavefunction in the Born-Oppenheimer approximation, but Eq. (1) is exact, as proven in the earlier works, and the equations that the two parts $\chi$ and $\Phi_{R}$ satisfy differ from that in the Born-Oppenheimer equation. These equations may be derived by inserting Eq. (1) into the Dirac-Frenkel action and applying the variational principle Abedi et al. (2010, 2012, 2013b); Alonso et al. (2013), and, for one nuclear degree of freedom and one electronic degree of freedom have the form

[TABLE]

[TABLE]

where time-dependent potential energy surface (TDPES) $\epsilon(R,t)$, vector potential $A(R,t)$, and coupling operator $\hat{U}_{\rm en}$ are given by

[TABLE]

with $\langle...\rangle_{r}$ indicating an inner product only over the electronic space. The potentials $V_{\rm ext}^{e}$ and $V_{\rm ext}^{n}$ represent potentials that are externally applied to the coupled system of electrons and nuclei, e.g. a laser field, and $\hat{H}_{\rm BO}$ is the usual Born-Oppenheimer Hamiltonian consisting of electronic kinetic energy, electron-electron, electron-nuclear, and nuclear-nuclear Coulomb potentials. Note that atomic units are used throughout, and $M$ denotes the effective mass of the nuclear coordinate. Under the phase-transformation of Eq. (II), the potentials $\epsilon$ and $A$ transform in a gauge-like way:

[TABLE]

The TDPES can be separated into a gauge-independent (GI) term and a gauge-dependent (GD) term, $\epsilon(R,t)=\epsilon_{\rm GD}(R,t)+\epsilon_{\rm GI}(R,t)$ where $\epsilon_{\rm GD}(R,t)=\langle\phi_{R}(t)|-i\partial_{t}\phi_{R}(t)\rangle$ and $\epsilon_{\rm GI}(R,t)$ are the remaining terms of Eq. (6). The three terms in Eqs. (6)– (8) capture the entire coupling of the electrons and nuclei; through them, the solution of the electronic equation depends on the nuclear wavefunction, and the solution of the nuclear wavefunction depends on the electronic wavefunction. Refs. Abedi et al. (2010, 2012, 2013a); Agostini et al. (2015); Min et al. (2014); Requist et al. (2016); Curchod and Agostini (2017) (for example), have emphasized the significance of $\epsilon$ and $A$ in yielding the exact nuclear dynamics; in contrast to the Born-Huang expansion, which is also exact, here the exact nuclear dynamics is achieved with a single potential energy surface rather than an infinite number of them, and a single vector potential coupling rather than an infinite number of first and second-order couplings. Wavepacket branching and decoherence are captured, where a central role is played by $\hat{U}_{\rm en}$.

III A Priori Numerical Considerations

Here we make some observations that will inform our attempt to solve Eqs. (4) – (8) self-consistently for a model system.

III.1 Reality of potentials $\epsilon(R,t)$ and $A(R,t)$

The equation for the nuclear wavefunction is of Schrödinger-form provided the potentials $A(R,t)$ and $\epsilon(R,t)$ are real. The reality of the potentials ensures that the Hamiltonian in the nuclear equation is Hermitian and consequently that the norm of the nuclear wavefunction is preserved during the time-evolution. However, the potentials $A(R,t)$ and $\epsilon(R,t)$ are determined by the solution of the electronic equation, which not only does not have a Schrödinger form, but has a Hamiltonian to which notions of hermiticity do not technically apply (see also Sec. III.2). Nevertheless, norm-conservation in the electronic equation is still a meaningful concept; indeed it is particularly important for the electronic equation, given the role of the PNC of Eq. (2) in ensuring that the factorization is unique. Propagation under the exact EF equations does conserve the partial norm Alonso et al. (2013); Abedi et al. (2013b) and in this subsection we show how this directly ensures the scalar and vector potentials that drive the nuclear equation are real. When the equations are discretized in a numerical simulation, additional considerations need to be made.

To see why the partial norm condition of the electronic wavefunction assures us that the scalar and vector potentials in the Schrödinger equation for the nuclear wavefunction are real, we first observe that

[TABLE]

where $\lambda$ could represent $R$ or $t$, i.e. a parameter of the electronic wavefunction. So the expectation value on the right-hand-side will be purely real if the PNC is satisfied. Taking $\lambda=R$, the expectation value appearing on the right-hand-side becomes the vector potential $A(R,t)$ and hence PNC condition ensures $A(R,t)$ is real. Taking $\lambda=t$, the expectation value on right-hand-side becomes the gauge-dependent part of the scalar potential, $\epsilon_{\rm GD}(R,t)$, hence the PNC condition ensures $\epsilon_{\rm GD}(R,t)$ is real. The first two terms in the gauge-independent part $\epsilon_{\rm GI}(R,t)=\langle\phi_{R}|\hat{H}_{\mathrm{BO}}+\hat{V}^{e}_{\mathrm{ext}}+\hat{U}_{\rm en}|\phi_{R}\rangle_{r}$ are real due to the operators $\hat{H}_{\mathrm{BO}}$ and $\hat{V}^{e}_{\mathrm{ext}}$ being Hermitian, while the last term can be shown to be $\langle\phi_{R}|\hat{U}_{\rm en}|\phi_{R}\rangle_{r}=\left(\langle\nabla_{R}\phi_{R}|\nabla_{R}\phi_{R}\rangle_{r}-A^{2}\right)/2M$ upon making use of the definition of the vector potential to cancel some terms. Thus, we have shown that the potentials $\epsilon(R,t)$ and $A(R,t)$ appearing in the Schrödinger equation for the nuclear wavefunction are real, provided that the PNC is satisfied.

However, while the above arguments hold for the exact (continuous) case, there is no guarantee they will always hold for a numerical solution because of discretization. For a given discretization, the analog of the relation Eq. (10) requires a certain additional choice of the discrete points at which to evaluate the integral. For example, consider the simplest discrete version of Eq. (10) taking $\lambda=R$: letting $\phi_{iJ}^{(n)}$ represent the conditional electronic wavefunction with $i,J,n$ as the index for $r,R,t$ respectively,

[TABLE]

The term inside the parentheses in the last line corresponds to $\langle\phi_{R}(r,t)|-i\nabla_{R}\phi_{R}(r,t)\rangle_{r}=A(R,t)$. As in the continuous case the fact that the norm is independent of $R$ ensures that $A(R,t)$ is real, but here only provided one takes the midpoint in $R$-space to evaluate $\langle\phi_{R}(r,t)|$, i.e. the average $(\phi^{(n)*}_{i,J+1}+\phi^{(n)*}_{i,J})/2$.

One can do the same calculation with time indices:

[TABLE]

Here the term in parenthesis is $\langle\phi_{R}(r,t)|-i\partial_{t}\phi_{R}(r,t)\rangle_{r}=\epsilon_{\rm GD}$. Again, as in the continuous case, PNC implies $\operatorname{Im}(\epsilon_{\rm GD})=0$ provided $\langle\phi_{R}(r,t)|$ is evaluated precisely at the midpoint in time.

In practise, we used higher-order discretizations in both space and time (see shortly) for which a more complex additional condition would be required to guarantee that satisfaction of the PNC implies reality of the potentials. We found, however, that for the time-steps and spatial-grid we used that neglecting this correction did not yield a large error, i.e. that any imaginary part in the potentials was small. A cruder fix that we tried was to simply to redefine the potentials, e.g. $A(R,t)\to\mathrm{Re}A(R,t)\;{\rm and}\nabla_{R}A(R,t)\to\mathrm{Re}\nabla_{R}A(R,t)$, but this did not appear to make much difference in the calculations for which we were able to propagate for a significant amount of time.

III.2 Hermiticity and partial norm conservation

Having discussed the subtleties in ensuring the reality of the potentials in the numerics assuming we have a propagation scheme that conserves the PNC, we now turn to what such a propagation scheme should be.

As mentioned in Sec. III.1, although the equation for the nuclear wavefunction is of Schrödinger-form, for which efficient and and accurate propagation schemes have been well-studied, the equation for the electronic equation is not. The form of the coupling potential $U_{en}[\phi,\chi]$ is unprecedented in the literature: it is non-linear and also has operators which act on the parametric dependence of the conditional electronic wavefunction. The domain of this operator is outside the electronic Hilbert space associated with a fixed nuclear configuration; rather, it connects electronic Hilbert spaces associated with neighboring nuclear coordinates. As a result it is not possible to define a Hermitian conjugate, a feature that is only compounded by its non-linearity; the notion of operator Hermiticity is not meaningful for general non-linear operators while a notion of Hermiticity can be rescued for such operators however in the context of expectation values, as in Ref. Schwartz (1997).

One of the most salient consequences of Hermiticity of the Hamiltonian in the usual Schrödinger equation, is that it guarantees norm-conservation. As mentioned in Sec. III.1, non-linear Hamiltonians may still preserve the norm, and the Hamiltonian in the electronic equation does. This is particularly important for uniqueness of the EF approach, and for ensuring that the coupling potentials appearing in the Schrödinger equation for the nuclear wavefunction are real.

A commonly-used time propagator for evolving the Schrödinger equation is the Crank-Nicolson (CN) propagator Crank and Nicolson (1996), which conserves the norm when the Hamiltonian is Hermitian. Considering a general discretized wavefunction $\psi^{(n)}_{i}$, with $n,i$ being the time and spatial index respectively, CN yields for the projection of the time-evolution $\langle\psi(t)|-i\partial_{t}\psi(t)\rangle$,

[TABLE]

which is real for Hermitian Hamiltonians $H$ (often evaluated at the midpoint in time for time-dependent and/or non-linear $H$).

We now consider the application of CN for the electronic equation. We start by expanding Eq. (4) and grouping terms two by two in the following way:

[TABLE]

such that each line contains two terms that cancel each other if one takes the expectation value over the conditional wavefunction. Now, we take the inner product with $\langle\phi_{R}(r,t)|$ to obtain the analog of Eq. (13). Then for partial norm conservation one only has to show that $\operatorname{Im}(\epsilon_{GD})=0$, with the $\epsilon_{GD}$ defined in Eq. (12). But in this case it is more complicated than the usual case in Eq. (13) as the zero imaginary part does not come from Hermiticity of the matrix. Rather, it results from exact cancellation of terms when the equation is projected on $\langle\phi_{R}(r,t)|$. For example $\langle\phi_{R}(t)|-\frac{1}{M}\frac{\nabla_{R}\chi(R,t)}{\chi(R,t)}\cdot\nabla_{R}|\phi_{R}(t)\rangle_{r}$ should cancel exactly $\langle\phi_{R}(t)|i\left[\frac{1}{M}\frac{\nabla_{R}\chi(R,t)}{\chi(R,t)}\cdot A(R,t)\right]|\phi_{R}(t)\rangle_{r}$. This obviously cancels in the continuous case from the definition of $A(R,t)$ but in the discrete case we have to define each element carefully. For example, if one defines

[TABLE]

to conserve the norm in a Crank-Nicolson fashion then we have to use an alternative definition of $A$, that we call $\tilde{A}(R,t)$ here. We define $\tilde{A}$ as:

[TABLE]

with ${\cal{N}}(R,t)$ the partial norm of the time midpoint:

[TABLE]

which is not equal to one. Choosing this definition does not enforce that $\tilde{A}(R,t)$ is real though it is equivalent to the definition following Eq. (11) when both $\Delta t$ and $\Delta R$ go to zero. Using the definition of Eq. (16) and renormalizing the average values of the form $\langle\phi_{R}(t)|O|\phi_{R}(t)\rangle_{r}$ that appear in Eq. (14) by ${\cal{N}}(R,t)$ would allow $\operatorname{Im}(\epsilon_{GD})$ to be zero and the partial norm to be conserved during the propagation. The main problem with this approach is dependence of many quantities on the next time-step with $\phi^{(n+1)}_{i,J}$ which would create a highly non-linear equation to solve during the propagation.

The analysis above demonstrates the challenge of discretizing the EF equations in a way that remains faithful to the norm-conservation aspects of the exact equations. In an alternative approach, one could consider instead discretizing first, and then factorizing, however this presents it own challenges in identifying separate equations with clear definitions for the potentials; some preliminary work in this direction is presented in Appendix A.

In our numerical investigations, we use a standard RK45 scheme (see Sec. III.5), with a small enough time-step and grid spacing that errors in partial norm conservation were small over the duration of our simulation. We found that enforcing a crude “norm-conservation correction” in our propagation in the form of

[TABLE]

generally improved the stability and accuracy of the dynamics.

III.3 EF expansion in the Born-Oppenheimer basis

Generally, a numerical solution of Eqs. (4) – (8) requires spatial grids for the nuclear coordinate $R$ and for the electronic coordinate $r$. Since electronic wavefunctions vary on the scale of tenths of an Angstrom while extending over a few Angstroms at least (more when the system evolves to more highly excited states or partially ionizes), the size of the electronic grid can limit the efficiency of the time-propagation. In some cases, the numerics can be made more tractable by using a basis for the electronic equation; the choice of an optimal basis depends on the physics of the problem, but for cases where external perturbations are weak and the system is initially well-described by one or a few BO states, expanding the conditional electronic wavefunction in a truncated BO basis can make the numerical simulation simpler and faster.

The equation of motion for the electronic system then turns into coupled equations for basis coefficients, $C_{j}(R,t)$, where

[TABLE]

with the adiabatic BO states satisfying $\hat{H}_{BO}(R,r)\phi^{j,BO}_{R}(r)=\epsilon^{j}(R)\phi^{j,BO}_{R}(r)$, with $\epsilon^{j}(R)$ the $j$th BO potential energy surface. Inserting this into the equation of motion for the conditional electronic wavefunction leads to (see Appendix (B) for details)

[TABLE]

where: the $i$th-projected $U_{\rm en}$ is

[TABLE]

the non-adiabatic coupling terms are

[TABLE]

and the remaining potential terms are

[TABLE]

The time evolution of the full molecular system then is equivalent to self-consistently evolving Eq. (20) with the equation for the nuclear wavefunction, Eq. (5).

III.4 Choice of gauge

While the EF equations are invariant under the gauge transform Eqn. (II), the stability properties may not be. It is possible that a ‘best gauge’ exists for the propagation of these equations, but for now we simply choose the gauge such that

[TABLE]

i.e., the gauge-dependent part of the time-dependent potential energy surface is set to zero.

The gauge condition Eq. (III.4) is not guaranteed to hold in a numerical solution and we do not enforce it explicitly. We could attempt to impose the gauge condition in a first-order way, replacing

[TABLE]

However we found that this does not appreciably extend the integration time with the propagation method we used ($<5$a.u. difference), even if it does ensure the system remains in the correct gauge for longer.

III.5 Propagation scheme

An important choice for numerical time-evolution is the propagation algorithm. Much depends on this choice and, for general non-linear systems, it is not always clear why one might be more suitable than another. Since the equation for the nuclear wavefunction has the time-dependent Schrödinger form, the Crank-Nicolson (CN) propagator is the preferred choice given its stability and norm-conserving properties. As discussed extensively in Sec. III.2, the electronic equation is not of this form and CN may not be the best choice due to the non-linearity. In our numerical investigations, we use two propagation schemes. For the majority of the presented results we evolve the electronic system with the fourth-order Runge-Kutta (RK45) scheme, in which the $R$ and $t$-dependent terms vary between the internal steps. Being a purely explicit integrator it is known to be only conditionally stable, and so to investigate RK45’s robustness to this problem we also test the simplest purely implicit scheme, the backwards Euler (BE) method, as well as with the mixed explicit-implicit CN scheme.

The inconsistency in the order of accuracy in this comparison ($O(\Delta t^{5})$ and $O(\Delta t^{2})$ for RK45 and BE/CN respectively) is somewhat intentional in that it allows us to investigate which is more important for a stable norm-conserving solution: accuracy or an implicit solver. For example, an implicit solver may have absolute stability over a large range of parameters but a low order method may expedite its exit from this region, whereas a more accurate method may be able to remain in its island of stability longer. A direct but preliminary analysis of the stability of the EF equations is carried out in Sec. V.

For the two implicit propagation schemes, BE and CN, a matrix solver is required: both BE and CN can be cast into the form ${\cal{M}}\cdot v=b$ for solution vector $v$. To avoid constructing matrices in memory, even when sparse, this equation is solved iteratively using the biconjugate gradient stabilized method (BiCGSTAB) van der Vorst (1992). This method assumes that $\cal{M}$ is not a function of $v$, and so requires linearization of the underlying equation. This is done by not updating the functions $A(R,t)$ and $\epsilon(R,t)$ (which are non-linear in the electronic coefficients) between steps of the BiCGSTAB procedure, leaving only the coefficients themselves and their gradients, all linear functions, free to vary.

III.6 Spatial discretization

The derivative operators, $\nabla_{R}$ and $\nabla_{R}^{2}$ are represented by five-point central finite-difference stencils, where the corresponding one-sided stencils are used at the boundaries. To prevent oscillations appearing at the tails of the nuclear wavefunction, $\chi$ and its first and second derivatives are set to zero at the boundaries.

IV Numerical Results

To investigate the propagation of the EF equations, we choose the one-dimensional Shin-Metiu system Shin and Metiu (1995); Abedi et al. (2013a); Agostini et al. (2015), which is a model of proton-coupled electron-transfer with one nuclear and one electronic coordinate. This model system has been instructive for studying different methods for non-adiabatic dynamics; in particular, the potentials that arise in the EF formalism have been studied and analyzed in this model, uncovering universal features such as steps in the TDPES bridging BO surfaces after a non-adiabatic event Abedi et al. (2013a); Agostini et al. (2015). In that work, the full TDSE was solved exactly, finding $\Psi(r,R,t)$, and extracting the vector and scalar potentials $A(R,t)$ and $\epsilon(r,t)$ via inversions. In contrast, here, we endeavor to directly solve the EF equations, with the potentials emerging on the fly during the propagation.

The Shin-Metiu model Hamiltonian corresponds to an electron and ion moving between two fixed ions a distance $L=19.0a_{0}$ apart:

[TABLE]

where the parameters, $M=1836$ a.u., $R_{f}=5.0a_{0}$, $R_{l}=3.1a_{0}$, $R_{r}=4.0a_{0}$. Fig. 1 shows the lowest two BO surfaces and the corresponding non-zero non-adiabatic couplings as defined in Eq. (22). The softened interparticle interactions in this model avoids any possible problems with Coulomb singularities Jecko et al. (2015).

The higher BO surfaces have a large enough energy separation from the first two such that, beginning with occupation only on the lowest two surfaces leads to only these two surfaces ever being occupied, at least for the duration of time that we consider. So we take a truncated basis of these two adiabatic levels for the electronic system. We take the same initial conditions as used in the earlier work Abedi et al. (2013a); Agostini et al. (2015), with the system starting on the upper BO electronic surface, $C_{1}(R,0)=0$, $C_{2}(R,0)=1$ $\forall R$, and with a gaussian nuclear wave function defined by $\chi(R,0)=\exp[-(R-R_{0})^{2}/2\sigma^{2}]$ with $R_{0}=4$ and $\sigma^{2}=1/5.7$.

Taking into account the considerations of Sec. III, we now embark upon the exact self-consistent solution of the EF equations. A grid of $2000$ $R$-points is used, corresponding to $\Delta R=0.009$, with a time step of, unless otherwise stated, $\Delta t=0.1$ a.u. We will compare the results of the EF simulation with the exact solution of the full two-dimensional TDSE with the CN propagator and the same finite difference stencils as above; here we use the same time-step and $R$-grid, and a grid-spacing of $0.3$au for the electronic coordinate.

A straightforward implementation of the CN propagation for $\chi$ together with the RK45 scheme for $C_{i}(R,t)$ (including the norm-conservation correction Eq. (18)), unfortunately becomes unstable and unbounded to the point of failure after a very short time, $t_{\mathrm{max}}=5$a.u., short enough that negligible dynamics have occurred. When we turned all the coupling terms to zero so that the equations reduced to BO dynamics, the instability vanished, and accurate propagation was obtained for long times, giving identical nuclear and conditional electronic densities as those obtained from the full solution of the electron-nuclear TDSE with zero coupling. This verifies that the problem lies in the coupling terms, and, in particular, the term $\nabla\chi/\chi$ appears to be problematic.

In the following, we describe different avenues to “tame” the instabilities, with varying degrees of success. Considering the best case scenario for different propagators, the maximum times for the simulation before it exploded, were as follows: 370 au for BE, 480 au for CN, and 1490 au for RK45. Below, the results for the RK45 simulation are studied.

IV.1 Masking the coupling terms

The instabilities are related to oscillations that initially develop away from the region of appreciable nuclear density before propagating inwards. In these regions, since $Re\nabla\chi/\chi=\nabla|\chi|/|\chi|$ can become unreliable, and, given that when the nuclear density is very small, the conditional electronic wavefunction has limited physical significance anyway, we apply a mask function which smoothly sets the coupling terms to zero far from the physical region.

The function used to generate the mask is defined in Appendix C, and consists of a smooth step on either side of the physical region. The mask ‘tracks’ the density throughout the simulation: we define the center of these steps to be the points at which the nuclear density drops below some threshold $\kappa$ as we approach the physical region from the left and from the right. An example is depicted in Fig. (2) along with the density. Unless otherwise specified, the values of $\kappa$ and $w$ indicated in the caption are used.

Outlined in the following subsections, are different choices as to which terms to apply the mask and the results on the dynamics.

IV.1.1 Masking $U_{\rm en}$($t_{\mathrm{max}}=610$a.u. with $\kappa=10^{-10}$)

Here we apply a mask to the entire coupling term $U_{\rm en}(r,R)$, so that, for all $r$, at the values of $R$ that have appreciable nuclear density (see above) the conditional electronic equation is as in Eq. (4) while book-ending either side of this region with purely Born-Oppenheimer dynamics for the conditional electronic wavefunction. The nuclear equation Eq. (5) is left unaltered.

With the mask parameters of Fig. 2, the calculation fails already at $t_{\mathrm{max}}=610$a.u. Increasing the value of $\kappa$ allows us to propagate a little further in time, and Figure 3 shows the result of propagation with this masked $U_{\rm en}$, depicting the nuclear density and TDPES at $t=790$au for two values of the mask edge tolerance $\kappa$. This shows that despite the best efforts of the mask on $U_{\rm en}$ to keep the influence of noise in the coupling terms in the asymptotic regions at bay, errors still propagate into the region of non-trivial density causing the simulation to fail. The level of the rapid oscillations in $\nabla\chi/\chi$ that reaches the area of significant density is apparently not strongly affected by the width of the mask, even for the constrictive case of $\kappa=10^{-3}$. A mask on $U_{\rm en}$ alone is is not sufficient to prevent errors from accumulating in the simulation, and any such errors are not confined to the periphery.

IV.1.2 Masking $\partial_{t}C_{i}(R,t)$: ($t_{\mathrm{max}}=340$a.u. with $\kappa=10^{-10}$)

An alternative way to mask the entire coupling term in the electronic equation is to directly mask the time-derivative of the electronic coefficients, that is, mask the entire right hand side of Eq. (20). Figure 4 shows the result of masking the $\partial_{t}C_{i}(R,t)$ functions given in Eqn. (21). Fig. 4 shows that this when the mask is placed on $\partial_{t}C_{i}(R,t)$ oscillations and spikes exist well inside the masked region even at relatively early times, and the simulation fails before appreciable density has reached the avoided-crossing region.

IV.1.3 Masking $\nabla\chi/\chi$ and $A(R,t)$ : ($t_{\mathrm{max}}=1490$a.u. with $\kappa=10^{-10}$)

An alternative to the above two approaches is to mask $\nabla\chi/\chi$ directly, as it is known to be a source of noise and instability early on. In fact, applying the same masking procedure as done above greatly improves the stability of the simulation and allows propagation long enough, to a time of $t_{max}=1220$ a.u., to witness the non-adiabatic event of wavepacket splitting. Additionally, we found that this simulation could be extended further by placing a mask on the vector potential $A(R,t)$ and its divergence $\nabla\cdot A$ at fixed points near boundaries We use the same mask function as is placed on $\nabla\chi/\chi$, however in this case we fix the start and the end of the mask to be a fixed distance, 100 grid points or 0.9 a.u., from the boundaries. This additional mask on the vector potential did not make a significant difference in the other cases above, but here, masking both $A,\nabla\cdot A$ and $\nabla\chi/\chi$ allowed us to propagate to $t_{max}=1490$ a.u. We therefore define a moving mask on $\nabla\chi/\chi$ with $\kappa=10^{-10}$ and $w=50$, as well as a fixed mask on the edges of the $A(R,t)$ and $\nabla\cdot A(R,t)$ as the best case scenario which we will refer to repeatedly in what follows, and what was plotted in Fig. 3 in black.

IV.1.4 Discussion on the best case scenario

This solution is depicted in Fig. 5 at several times, showing the exact TDSE and EF nuclear densities, the two adiabatic surfaces, the TDPES, and the mask. In all panels the solution agrees well with the direct TDSE solution. One sees the diabatic nature of the TDPES as the density passes through the avoided crossing region, followed by bridging of the two BO surfaces after passage through the avoided crossing region Abedi et al. (2013a) that accompanies the splitting of the nuclear wavepacket. We note that in the chosen gauge the gauge-dependent part of the TDPES is zero, so no piecewise off-set is seen Agostini et al. (2015).

Noise in the TDPES becomes clearly visible from 1000a.u. onwards, and appears to arise initially where the density is small near the mask boundary; this is where $Re\nabla\chi/\chi=\nabla|\chi|/|\chi|$ is largest. We found this also with different values of $\kappa$ which broadens or narrows the mask.

To investigate whether simply having large values of $\nabla\chi/\chi$ instigates noise we inserted an analytic form for $\nabla\chi/\chi$ into the electronic equation to see if unstable behavior persists. Naturally, if the term we use bears no relation to the true $\nabla\chi/\chi$ then we are no longer solving the physical problem and any observations we make may have limited scope and utility. Still, it is instructive to understand the numerical properties of this system: with no mask, if $\nabla\chi/\chi$ is replaced by an analytic form, does the solution of the coupled EF equations still fail quickly? We fix $\nabla\chi/\chi$ to be purely real and have some artificial form. With this we evolve the electronic system, which is now totally decoupled from the nuclear system, and find that the simulation gives a catastrophically noisy TDPES almost immediately ($<10$au). Fig. 6 shows simulations, with no masks, for three different forms of $\nabla\chi/\chi$ fed in to the electronic equation: (i) $\nabla\chi(R,t=0)/\chi(R,t=0)$, (ii) $\nabla\chi(R,t=0)/\chi(R,t=0)+C$ so that it has the same slope but is vertically shifted with $C$ such that the maximum magnitude is on the right-hand boundary, and (iii) a Gaussian with a maximum at the peak density, which has the least physical resemblance to the actual function. From the results, we immediately observe that one can induce instability in the electronic system at a particular $R$ simply by having $\nabla\chi/\chi$ have a large value there. Thus this numerical exploration strongly suggests that a key source of instability in propagation of the exact EF equations is large $\nabla\chi/\chi$. We come back to this point in the preliminary mathematical analysis of the equations in Sec. V.

IV.2 Reducing the time step

A common cause of noise and instability in numerical simulation is violation of a Courant-Friedrichs-Lewy (CFL) condition which provides an upper bound on the time step for a given grid spacing Courant et al. (1967). To investigate this we reduced the step size by an order of magnitude for the “best case scenario” simulation with the masked $\nabla\chi/\chi$ and $A$. Surprisingly, this simulation failed earlier, at 1170a.u. The TDPES and Re$(\nabla\chi/\chi)$ at this time are plotted for both time-steps in Fig. 7, showing that indeed the smaller $\Delta t$ choice develops oscillations and sharp features where the larger $\Delta t$ choice has none. A more detailed analysis of this situation is given in Sec. (V).

V Stability analysis

The fact that reducing the time-step causes the numerical solution to become unstable at shorter times indicates that the stability properties of the system of equations, Eqs. (4) – (8), are unusual and highly non-trivial. The electronic equation has an unprecedented form, whose mathematical and numerical properties are unknown. Do numerical methods exist in which the solution is stable and converged to the true solution? Here we only begin to scratch the surface of the question of whether, once discretized in time and space, the equations are numerically stable, with standard propagation methods. We do not address here the separate question of whether the exact EF equations are mathematically well-posed.

The nuclear equation, Eq (5), is a time-dependent Schrödinger equation for which stability properties are well-known for potentials that are possibly time-dependent but pre-determined (i.e. not determined by a self-consistent coupled solution as in the EF case). Specifically, we consider the equations for the electronic coefficients, $C_{i}(R,t)$, which we place into a vector, $(C_{1}(R,t),C_{2}(R,t)...C_{k}(R,t))$ ($k$ truncated to 2 in our model system of Sec IV but now we keep it general). Then we write Eq. (20) as a matrix equation, as if it was a linear equation,

[TABLE]

where underline denotes a matrix, and we have separated out the terms as follows:

$\underline{f}_{0}$ has a Schrödinger structure, i.e. matrix proportional to $\delta_{ij}$ with elements $\epsilon^{i}({\bf R},t)-\epsilon({\bf R},t)-\nabla^{2}/2M$. This is diagonal in sense of not coupling different $C_{i}$, but it is not diagonal once $C_{i}$ is resolved in $R$ due to the Laplacian, i.e. it would be block-diagonal when writing each coefficient out on a grid $(C_{1}(R_{1},t),C_{1}(R_{2},t)...C_{1}(R_{N},t),C_{2}(R_{1},t),C_{2}(R_{2},t).....)$

$\underline{f_{1}}$ contains all the other purely multiplicative terms (including those coupling the different $C_{i}$’s via $d_{ij}$ and $d^{2}_{ij}$).

$\underline{g}$ contains terms involving the first derivative $\nabla_{R}$ and is not diagonal in $C_{i}$.

The actual equations for the coefficients are however not linear since $A$, $\epsilon$ depend on $C_{i}$ (see Eq. (23)), but to simplify the analysis for a preliminary examination of stability, we consider these as just some functions of $(R,t)$, and further, that they are determined ahead of time, i.e. they are not produced by some iteration with the nuclear equation, as they would be in actuality. With this (big) simplification, what can we say about stability of numerical solutions of this system?

Discretizing in time, we write Eqs. (V) as

[TABLE]

where $n=0,1,2..$ labels the time-step. The system will be numerically stable if the magnitude of the eigenvalues of the amplification matrix G are bounded from below by 1. We will consider the simplest schemes: forward and backward Euler. From experience with the TDSE and the $\underline{f}_{0}$ “backbone” of our equation, we do not expect that forward Euler will be stable, but it is nevertheless instructive to recall why.

The forward Euler (FE) method is defined by

[TABLE]

and so

[TABLE]

If we take just the TDSE-type terms, $\underline{f}_{0}$, this is unconditionally unstable for any finite time-step, because then $\underline{K}=\underline{f}_{0}$, being Hermitian, has real eigenvalues, which leads to the (complex) eigenvalues of $\underline{G}_{\mathrm{FE}}$, $\lambda_{\rm FE}=1-i\Delta t\lambda_{K}$ where $\lambda_{K}$ are the (real) eigenvalues of $\underline{K}$, having magnitude $|\lambda_{\rm FE}|=\sqrt{1+\Delta t^{2}\lambda_{K}^{2}}$, greater than 1. With all the terms included in $\underline{K}$, it will also not be stable in general, especially given that at the early times the dynamics occurs on the first excited BO surface, with $\underline{f}_{0}$ dominating $\underline{K}$, until appreciable density approaches the avoided crossing.

We next consider an implicit method, as is usually done for TDSE, the simplest of which is the backward Euler (BE) method:

[TABLE]

Noting that eigenvalues of the matrix inverse are the inverse of the eigenvalues of the matrix, we see that for $\underline{K}=\underline{f}_{0}$, the eigenvalues have magnitude $|\lambda_{\rm 0,BE}|=\frac{1}{\sqrt{1+\Delta t^{2}\lambda_{K}^{2}}}\approx 1-\Delta t^{2}\lambda_{K}^{2}/2$ which is less than one for any non-zero $\Delta t$, so this is unconditionally stable, as is well known for TDSE-form equations. For the remainder of this section, we consider only BE.

Now we ask how the situation changes when the coupling terms in $\underline{f}_{1}$ and $\underline{g}$ are included. Because these terms are not Hermitian, the eigenvalues of $\underline{K}$ are no longer guaranteed to be real. Their imaginary part, once multiplied by $i\Delta t$, yields a correction to the real part of the eigenvalue of $G_{\rm BE}$, decreasing or increasing it away from 1, as well as contributing to the imaginary part. To make the analysis tractable, we adopt two simplifications. First, we truncate the number of coefficients to two (as was relevant in the numerical example of the previous section). Second, we use a minimal spatial grid for the finite-difference derivatives, taking a three-point centered stencil for the Laplacian and a two-point centered stencil for the first derivative. Third, we consider all the other $R$-dependent functions ($\epsilon$, $\epsilon_{i}$,$d_{ij}$, $d^{2}_{ij}$) as spatially-uniform; which is not a bad approximation within the 3$\times$3 spatial-grid truncation if our $R$-grid is fine enough. Although these approximations may seem drastic, they are valuable in that they straightforwardly lead to a preliminary analysis of the stability properties of the electronic equation. With these three approximations in hand, we write

[TABLE]

where $A_{i}$ and $B_{i,j}$ are tridiagonal matrices in the basis $\left[n-1,n,n+1\right]$ with elements

[TABLE]

where

[TABLE]

and

[TABLE]

For small $\Delta t$, we find only one non-trivial eigenvalue, which in this small-$\Delta t$ limit, is

[TABLE]

If this is greater than 1 for some $(R,t)$ then the solution locally is exponentially increasing, leading to an instability; a similar result holds for the Forward Euler. The local stability, for small time-steps $\Delta t$, depends then in a subtle way on the sign and size of

[TABLE]

and if we were to try to define a CFL-condition, it would be a spatially- and time- dependent condition, depending on the self-consistent solution of $\chi$ and $\phi$ at each time-step and spatial $R$-point. This result relates back to our numerical observation that instability appears to be triggered when $\nabla|\chi(R,t)|/|\chi(R,t)|$ gets large in magnitude. If $\Delta t$ is chosen small enough, it is necessary at least for the term in Eq. (37) to be negative, but there is no obvious physical reason why these terms should be negative everywhere.

The analysis above relied heavily on simplifications, exploring only the simplest finite-difference and time-propagation schemes. It is quite possible that a more sophisticated numerical scheme can be developed that is unconditionally stable. The results of the analysis should therefore not be a deterrent to finding such a scheme, but rather serve to point out that the stability conditions are far more complex than in most time-propagation schemes we have encountered.

VI Conclusions and Outlook

The potentials and coupling operator in the EF equations have already provided much insight into fundamental aspects of coupled dynamics of subsystems, while a mixed quantum-classical scheme based on them has already proven to be a soundly-based approach for practical problems. The investigations here on the self-consistent numerical solution of the equations have fundamental interest but also a practical interest: it would be very useful in developing mixed quantum-classical approximations to test the effect of approximations made to the terms prior to making any classical or semiclassical approximation. However, this work has shown such an endeavor is very challenging. Care must be taken to ensure that features as fundamental as having real potentials and partial norm conservation are generally robust in the numerical scheme, yet the structure of the electronic equation shows the usual propagation methods will not respect this. These methods applied to propagate the EF equations become unstable rather quickly. We showed that the numerical instability can be tamed by using mask functions on the coupling term $\nabla\chi/\chi$, and the vector potential $A$, such that numerically-exact propagation of the EF equations in a model system can be achieved for a time long enough that non-adiabatic branching phenomena and wavepacket splitting can be observed. But even with the masks, instabilities kill the calculation soon after. The purpose of this work was not to present an exhaustive exploration of different algorithms but rather to illustrate some of the unforeseen challenges when attempting a direct solution. One cannot pinpoint one term that is the culprit causing the instability, as we saw instabilities remain when $\nabla\chi/\chi$ was replaced with a smooth, analytic, function. That there is a complex interplay of terms is also suggested by our preliminary stability analysis exploring the CFL-condition for an implicit propagation scheme. We hope that this paper will motivate further work to take these developments further for a fuller understanding of the stability of the EF equations and approximations.

Acknowledgements.

We thank Eric Cances for helpful discussions. We are grateful to the U.S. National Science Foundation grant CHE-1566197 (L. L. and G. G.), Department of Energy Office of Basic Energy Sciences, Division of Chemical Sciences, Geosciences and Biosciences grant DE-SC0015344 (N. T. M. and G. G. ), and the Research Corporation for Science Advancement Cottrell Scholar Seed Award (G.G.) for support of this work.

Appendix A Factorization Following Discretization

We saw in Sec. III.2 that discretization of the EF equations raises challenging issues for conventional time-propagators. Instead of factorizing and then discretizing, one can choose a certain discrete propagator for the full wavefunction first, and then factorize the discrete equation. The two equations obtained this way should be numerically equivalent to the first equation.

Here we start with the equation for $\Psi^{(n+1)}_{i,j}$ in the Crank-Nicolson scheme:

[TABLE]

We define the following compact notations:

[TABLE]

Now we write $\Psi^{(n)}_{i,j}=\chi^{(n)}_{j}\phi^{(n)}_{i,j}$ and expand the operators $\partial_{t}$ and $\hat{T}_{n}$ for the discretized wavefunctions. We write here the first and second order derivatives of a product of two arbitrary functions $g$ and $f$:

[TABLE]

[TABLE]

which can be written in a more compact way as

[TABLE]

and

[TABLE]

We see the notion of midpoints (Sec. III.1) appear again.

Using these equations we rewrite Eq. (38):

[TABLE]

where $\phi^{(n+1/2)}=\frac{\phi^{(n+1)}_{i,j}+\phi^{(n)}_{i,j}}{2}$. In the continuous case one would take the inner product with $\langle\phi_{R}(r,t)|$ to obtain the equation of motion for $\chi(R,t)$. In this case, the equivalent would be to multply Eq. (44) by $\phi^{(n+1/2)}/\cal{N}$ and integrate (sum) over $i$. Contrary to the coninuous case, in the discrete case the right-hand-side does not simplify and identifying quantities like $A(R,t)$ or rewriting in a familiar form with a $\hat{U}_{en}$ seems impossible.

Appendix B Equation of motion for $C_{j}(R,t)$

In the factorization $\Psi(r,R,t)=\phi_{R}(r,t)\chi(R,t)$ the equation of motion for $\phi_{R}(r,t)$ is given by

[TABLE]

with electronic hamiltonian given by

[TABLE]

In the following calculations for simplicity we neglect the external potential as it does not enter into the Shin-Metiu test case currently under consideration. Rearranging terms in Eq. 8, we write the coupling potential as

[TABLE]

where $A(R,t)$ is the vector potential given by $A(R,t)=\langle\phi_{R}(r,t)|-i\nabla_{R}\phi_{R}(r,t)\rangle$.

We now expand the electronic wavefunction in terms of Born-Oppenheimer (BO) basis vectors as

[TABLE]

where

[TABLE]

with $\epsilon^{j}(R)$ the $j$th BO surface. All spatial derivatives of $\phi_{R}(r,t)$ transform into spatial-derivatives of the basis, which, being time-independent, need only be computed once at the outset of the simulation, and are thus immune to error propagation over the course of the simulation.

For brevity we now drop the ‘BO’ superscript for the basis vectors and so henceforth all wavefunctions denoted by $\phi^{j}_{R}(r)$ are assumed to be BO vectors.

Projecting the electronic equation,

[TABLE]

onto the state $\langle\phi^{i}_{R}(r)|$ we obtain

[TABLE]

where we have used Eq. (49) and orthogonality of the basis vectors to collapse a number of sums.

To compute $\sum_{j}\langle\phi^{i}_{R}(r)|\hat{U}_{en}\left[C_{j}(R,t)|\phi^{j}_{R}(r)\rangle\right]$ we expand the derivatives of Eq. (47):

[TABLE]

Projecting onto state $\langle\phi^{i}_{R}(r)|$ we get

[TABLE]

Defining the non-adiabatic coupling vectors as in Eq. (22), we find Eq. (21).

Likewise, expanding the vector potential and time dependent potential energy surface in the BO basis yields Eqs. 23.

Appendix C Mask function

The mask function used in this work is defined by taking the product of two step functions, one stepping ‘up’ on the left hand side and the other ‘down’ on the right hand side. Defining $R_{m}$ to be the position of a single one of these steps, $w$ a measure of the steep-ness of the step (explain), and $s=1$ for a step up and $s=-1$ for a step down, the single-step function is defined by

[TABLE]

where

[TABLE]

Thus to mask away the edges of a function as done for $U_{en}$ and $\nabla\chi/\chi$ one multiplies those functions by $m(R,R_{L},1,w)\times m(R,R_{R},-1,w)$ where $R_{L}$ and $R_{R}$ are the left and right hand edges of the mask respectively.