An extension of the fewest switches surface hopping algorithm to complex Hamiltonians and photophysics in magnetic fields: Berry's phase and "magnetic" forces

TL;DR

This paper extends the fewest switches surface hopping algorithm to complex Hamiltonians to better model photoexcited molecules in magnetic fields, incorporating Berry's phase effects and magnetic forces.

Contribution

It introduces a novel extension of FSSH for complex Hamiltonians, including Berry's phase and magnetic forces, addressing limitations of the original algorithm.

Findings

Partial success in capturing Berry's phase effects

Incorporation of magnetic forces in surface hopping

Highlighting limitations of local FSSH for topological effects

Abstract

We present a preliminary extension of the fewest switches surface hopping (FSSH) algorithm to the case of complex Hamiltonians as appropriate for modeling the dynamics of photoexcited molecules in magnetic fields. We make ansatze for the direction of momentum rescaling and we account for Berry's phase effects through "magnetic" forces as applicable in the adiabatic limit. Because Berry's phase is a nonlocal, topological characteristic of a set of entangled potential energy surfaces, we find that Tully's local FSSH algorithm can only partially capture the correct physics.

Equations

$$\displaystyle H(\bm{r})=\left(\begin{array}[]{cc}V_{00}(\bm{r})&V_{01}(\bm{r})\\ V_{10}(\bm{r})&V_{11}(\bm{r})\end{array}\right)$$

$$\dot{\mathbf{r}}$$

$$\dot{\mathbf{p}}$$

$$g_{0\to 1}$$

$${\mathbf{F}}_j^{mag}$$

$${\mathbf{B}}_j$$

$${\mathbf{F}}_j^{mag}$$

$${\mathbf{n}}_{jk}$$

$$H$$

$$\theta$$

$$\varphi$$

$${\lambda }_0$$

$${\lambda }_1$$

$${\psi }_0$$

$${\psi }_1$$

$${\mathbf{d}}_{01}$$

$${\Psi }_0(\mathbf{r})（$$

$${\Psi }_1(\mathbf{r})（$$

$${\mathbf{F}}_1^{mag}$$

$${\mathbf{F}}_0^{mag}$$

$${\mathbf{n}}_{01}=\left(({\partial }_x\theta {)}^2{p}^x,({\partial }_y\varphi \sin \theta {)}^2{p}^y\right)$$

$$f(\eta )$$

$$\displaystyle\bm{n}_{01}=\left\{\begin{array}[]{ll}(1,0)&\mbox{ when }(\partial_{x}\theta)^{2}>(\sin{\theta}\partial_{y}\phi)^{2}\\ (0,1)&\mbox{ when }(\partial_{x}\theta)^{2}<(\sin{\theta}\partial_{y}\phi)^{2}\\ \end{array}\right.$$

$$p_1^y$$

$$p_0^y$$

$${\mathbf{p}}_1$$

$$=-i\hslash (\frac{C_1^{\ast }\nabla C_1}{C_1^{\ast }{C}_1}+⟨{\psi }_1|\nabla |{\psi }_1⟩)$$

$$f(\eta )$$

$$=\frac{1}{2}\left(||{\mathbf{d}}_R|{|}^2+||{\mathbf{d}}_I|{|}^2\right)+\frac{\cos 2\eta }{2}\left(||{\mathbf{d}}_R|{|}^2-||{\mathbf{d}}_I|{|}^2\right)-\sin 2\eta {\mathbf{d}}_R\cdot {\mathbf{d}}_I$$

$$\tan 2\eta =\frac{-2{\mathbf{d}}_R\cdot {\mathbf{d}}_I}{||{\mathbf{d}}_R|{|}^2-||{\mathbf{d}}_I|{|}^2}$$

$$\cos 2\eta \left(||{\mathbf{d}}_R|{|}^2-||{\mathbf{d}}_I|{|}^2\right)>0$$

Full text

I Introduction

Fewest switches surface hopping (FSSH)Tully (1990) has been a very powerful tool for simulating non-adiabatic dynamics over the last thirty years.Oberhofer et al. (2017); Habenicht and Prezhdo (2008); Nelson et al. (2011); Sterpone et al. (2009); Nelson et al. (2014) The basic idea of the FSSH algorithm is to run stochastic dynamics on electronic adiabats, with stochastic switches between adiabats to account for electronic relaxation; in the spirit of Pechukas’s force,Tully (1991); Pechukas (1969) one rescales momenta in the direction of the derivative coupling whenever a hop between surfaces occurs. The algorithm has been shown to successfully capture both the short time dynamics of non-adiabatic systemsLandry and Subotnik (2014); Jain and Subotnik (2015) as well as (their) long time equilibrium properties.Parandekar and Tully (2005) At the same time, the cost of FSSH is quite modestTully (1990); Jain et al. (2016); Barbatti et al. (2007). Of course, Tully’s algorithm has a few well-known shortcomings: $(i)$ the original algorithm did not treat wave packet separation correctly, and thus did not model decoherence;Schwartz et al. (1996); Prezhdo and Rossky (1997a); Wong and Rossky (2002a, b); Fang and Hammes-Schiffer (1999a, b); Hack and Truhlar (2001); Volobuev et al. (2000); Jasper and Truhlar (2005); Prezhdo and Rossky (1997b) $(ii)$ the algorithm does not treat recoherence correctly;Subotnik et al. (2016) $(iii)$ the algorithm does not include any nuclear quantum effects.Craig and Manolopoulos (2004); Jang and Voth (1999) Of the problems above, item $(i)$ has been discussed extensively in the literature and can largely be corrected; items $(ii)$ and $(iii)$ are largely intractable with classical, non-interacting trajectories.Donoso and Martens (2001); Donoso et al. (2003) Nevertheless, as a testament of the algorithm’s value, FSSH is routinely applied today to simulate non-adiabatic dynamics including photochemical processesMüller and Stock (1997); Landry and Subotnik (2014); Nelson et al. (2011), scatteringShenvi et al. (2009); Golibrzuch et al. (2014); Jain et al. (2016), and charge transfer in solutionFang and Hammes-Schiffer (1999a); Schwerdtfeger et al. (2014).

Interestingly, of all of the applications listed above, there is one glaring omission. To our knowledge, no one has yet used FSSH to study non-adiabatic dynamics for molecular systems with spin degrees of freedom in strong magnetic fields. More generally, to our knowledge, no one has yet extended the FSSH algorithm to treat complex (rather than real-valued) electronic Hamiltonians. When one considers such an extension, several obvious questions arise, including: $(a)$ how should one incorporate geometric phases semiclassically,Berry (1984); Baer (2006); Yarkony (1996a, b) given that geometric phase is a nonlocal, topological property? $(b)$ how should one choose the direction of momentum rescaling when the derivative coupling is complex, and there is no unique real vector to isolate? The answers are not obvious.

With this background in mind, the goal of this paper is to propose one possible set of answers and a possible extension of FSSH to the case of complex Hamiltonians. We will find that our current implementation of FSSH behaves reasonably well, though one clearly loses some accuracy when moving from the case of real to complex Hamiltonians. In particular, because of topological phase effects, we will show that obvious limitations arise for any algorithm (like FSSH) based on independent, spatially local and time local trajectories. This paper is structured as follows: In Sec. II, we introduce our several ansätze for the FSSH algorithm in the presence of a complex Hamiltonian. In Sec. III, we make clear our simulation details. In Sec. IV, we present our results. In Sec. V, we interpret our numerical results and give a simple explanation for how geometric phase effects appear in surface hopping, and we propose a general extension of the FSSH algorithm. Finally, in Sec. VI and Sec. VII, we summarize the paper and present some open questions, respectively. As far as notation is concerned, below we use bold characters (e.g. $\bm{r}$) to denote vectors, and we use plain characters (e.g. $H$) to denote either scalars or operators.

II Methods

II.1 Real Hamiltonians

Let us now briefly review the FSSH algorithm. As originally conceived, the FSSH approach is applicable to the case of real electronic Hamiltonians. Without loss of generality, consider a real two-by-two Hamiltonian (i.e. a Hamiltonian with two electronic states) of the form

[TABLE]

Here, $\bm{r}$ is a nuclear coordinate. To simulate semiclassical dynamics with quantum electronic states and classical nuclei, according to FSSH, one first diagonalizes the electronic Hamiltonian and computes adiabatic energies $E_{0}(\bm{r}),E_{1}(\bm{r})$, forces $\bm{F}_{0}(\bm{r}),\bm{F}_{1}(\bm{r})$, and derivative couplings $\bm{d}_{01}(\bm{r})$. Thereafter, one runs an ensemble of independent trajectories, initialized so as to correspond the correct Wigner distribution at time zeroWigner (1932); Subotnik et al. (2016); Landry et al. (2013); Subotnik et al. (2013); Kapral (2016). Each trajectory is evolved along a single adiabatic surface, with equations of motion:

[TABLE]

Here $j$ is the active surface for a given trajectory. Occasionally, trajectories switch from one surface to the other. For example, Tully proposedTully (1990) that a trajectory on surface [math] switches to surface $1$ with rate

[TABLE]

Here $\rho_{jk}\equiv c_{j}c_{k}^{*}$ are density matrix elements, and $(c_{0},c_{1})$ is the electronic wavefunction. Whenever a particle switches surfaces, in order to conserve energy, one rescales the momentum in the direction of the derivative coupling, $\bm{d}_{01}(\bm{r})$.Herman (1982) There are many existing references in the literature where one can learn more details of the FSSH algorithmTully (1990); Fang and Hammes-Schiffer (1999a); Barbatti (2011); Subotnik et al. (2016), beginning with Tully’s original paper. Tully (1990)

II.2 Complex Hamiltonians

At this point, we come to the heart of the matter. Consider a situation whereby a particle with spin interacts with a magnetic field and there are two possible electronic states. Because of the magnetic field, the electronic Hamiltonian will no longer be real-valued.Mead (1979); Yabushita et al. (1999); Domcke et al. (2004); Belcher (2011) Instead, the electronic Hamiltonian will be complex and therefore violate time reversibility. In other words, in Eq. 1, $V_{01}(\bm{r})$ can have both real and imaginary parts,Domcke et al. (2004) and $V_{10}(\bm{r})=V_{01}^{*}(\bm{r})$. For this situation, FSSH is not well defined and two obvious problems present themselves.

First, note that FSSH depends critically on the existence of adiabatic states. Now, it is well known that, in the presence of conical intersections, adiabatic electronic states cannot be globally defined, even for real electronic Hamiltonians.Domcke et al. (2004) Nevertheless, even though FSSH does not account for geometric phase, the algorithm is largely able to model dynamics through conical intersections, as has been documented in detail previouslyRyabinkin and Izmaylov (2013); Ryabinkin et al. (2014, 2017); Gherib et al. (2015). That being said, for the present case of a complex Hamiltonian, one must always worry: How should one best choose the sign of the wavefunctions, when the sign has a true complex phase and not just a plus/minus? And how should one best incorporate Berry’s phase effectsBerry (1984); Shankar (2012) semiclassically? 2. 2.

The second obvious question is: What is the (real-valued) direction for rescaling momentum? Obviously $Re(\bm{d}_{01})$ is not acceptable as this quantity depends on the choice of phase for the adiabatic electronic states. Furthermore, for a practical FSSH calculation, we must be able to compute this direction using only local information at a single nuclear geometry.

With these two questions in mind, we will propose a few simple and robust extensions of FSSH to complex Hamiltonians.

II.2.1 “Magnetic Force” Ansatz

As far as the changing (Berry) phase of the adiabatic electronic states, it is well known that the Berry curvature near the crossing region can be transformed into an effective magnetic field that is applicable in the adiabatic limit.Shankar (2012); Sakurai et al. (2014) Thus, to incorporate Berry’s phase effects into FSSH dynamics, we propose that, when a trajectory is moving on adiabatic surface $j$ near a crossing point, we will allow each FSSH trajectory to feel this extra “magnetic force”:

[TABLE]

Here, $\bm{B}_{j}$ is defined to be the Berry curvatureBerry (1984); Mead and Truhlar (1979)

[TABLE]

Substituting Eq. 5 into Eq. 4, and utilizing the identity $\bm{d}_{jk}=-\bm{d}_{kj}^{*}$, we find

[TABLE]

In the end, for an FSSH simulation moving along adiabat $j$, we will assume that the “magnetic” force $\bm{F}^{mag}_{j}$ should simply be added to the total adiabatic, Born-Oppenheimer force in Eq. 2. Note that $\bm{p}\cdot\bm{F}^{mag}_{j}=0$, so that this extra “magnetic” force does not break energy conservation, but rather turns the direction of momentum. Note further that this “magnetic” force disappears for the case of a real-valued Hamiltonian, where the derivative coupling $\bm{d}_{jk}$ is real. Interestingly, for a two state problem, Eq. 6 implies that $\bm{F}^{mag}_{0}=-\bm{F}^{mag}_{1}$.

II.2.2 Direction of Momentum Rescaling

In order to extend FSSH to the case of a complex Hamiltonian, we must find an appropriate direction for momentum rescaling, $\bm{n}_{jk}$, when a hop between adiabats $j\rightarrow k$ occurs. To be appropriate, this direction vector must satisfy at least three constraints: $(i)$ $\bm{n}_{jk}$ must be real; $(ii)$ $\bm{n}_{jk}$ should not depend on the phase of the derivative coupling $\bm{d}_{jk}$; $(iii)$ $\bm{n}_{jk}$ must reduce to $\bm{d}_{jk}$ when the complex part of the Hamiltonian is removed. Furthermore, we must be able to construct this direction with only local information at a single nuclear geometry; we cannot assume that we have any information about a global reaction coordinate.

With these constraints in mind, the following three ansätze for $\bm{n}_{jk}$ are possibilities:

• •

Method #1: “$Re(\bm{d}(\bm{v}\cdot\bm{d}))$”

Because the magnetic force is independent of phase, the following ansatz would appear reasonable:

[TABLE]

Note the strong connection between the magnetic force (Eq. 6) and $\bm{n}_{jk}$ here: According to Eq. 7, the real part of $\left[\bm{d}_{jk}\left(\frac{\bm{p}}{m}\cdot\bm{d}_{kj}\right)\right]$ would act as a direction for momentum rescaling while the imaginary part acts as a magnetic force that modifies motion along a given adiabat (see Eq. 6).

• •

Method #2: “$Re(e^{i\eta}\bm{d})$”

Another option for the rescaling direction $\bm{n}_{jk}$ is the real part of the derivative coupling with a robust phase factor. To this end, one can choose

[TABLE]

where for every coordinate $\bm{r}$, $\eta$ is chosen so as to maximize the vector norm $||Re(e^{i\eta}\bm{d}_{jk})||^{2}$. Note that, unlike Method #1, this ansatz for $\eta$ does not depend on any dynamical properties of a given trajectory.

• •

Method #3: “Average $\bm{d}$”

One last possibility is the averaged derivative coupling (divided by $2i$) 111Likely, $\bm{n}_{jk}=Re(\rho_{jk}\bm{d}_{kj})$ is another possibility.

[TABLE]

Like Method #1, this ansatz depends on the dynamics of a given trajectory. However, whereas Method #1 makes use of the nuclear momentum, Method #3 makes use of the electronic density matrix to construct the rescaling direction.

In practice, as shown in the Appendix A, Method #3 performs very poorly, 222Similarly, $\bm{n}_{jk}=Re(\rho_{jk}\bm{d}_{kj})$ does not perform well. and so below we will focus exclusively on Methods #1 and #2.

Throughout this paper, there is one nuance worth reporting. When running FSSH calculations, one needs to choose appropriate phases for eigenvectors. To choose these phases, one can use either $(i)$ eigenvectors computed on the fly, whereby the phase of a given set of eigenvectors are aligned with the eigenvectors at previous time step by “parallel transport” (i.e. $\left\langle\psi_{i}(t)\middle|\psi_{i}(t+dt)\right\rangle\approx 1$); or $(ii)$ analytical eigenvectors (see below in Eq. 12) for which a global phase is assigned (whenever possible). In our FSSH calculations, we find that as long as we initialize the system in a consistent fashion, we can use either phase convention, the difference between $(i)$ and $(ii)$ is negligible. For results below, all FSSH data are implemented using option $(i)$.

III Simulation Details

Consider a simple 2-D system with the following general Hamiltonian:

[TABLE]

For a simple model, we define the functions $\theta(x,y)$ and $\phi(x,y)$ to be:

[TABLE]

Here $A$, $B$, $W$ are constants. In Fig. 1, we plot the diabats, adiabats and derivative couplings. Note that the adiabats are completely flat which will make all FSSH results easier to interpret.

For this Hamiltonian, one may solve for the eigenvalues, eigenvectors, and the derivative couplings analytically:

[TABLE]

Note that, with the choice of adiabats in Eq. 12, $\bm{d}_{01}$ is composed of two components: a real component in the direction of the crossing ($\bm{\nabla}\theta$) and an imaginary component in the direction of the gradient of the phase of the diabatic coupling ($\bm{\nabla}\phi$). For all dynamics reported below, we initialize Gaussian wave packets on the upper surface

[TABLE]

Here $\bm{p}_{init}$ and $\bm{r}_{init}$ are the initial momentum and position, respectively; $\sigma$ is the spread of the initial wave packet over real space. For exact quantum calculations, the wave packets are propagated with the Schrödinger equation using the fast Fourier transform techniqueKosloff and Kosloff (1983). For the surface hopping algorithm, $10^{7}$ trajectories are sampled from the Wigner distribution corresponding to Eq. 13. Each semiclassical trajectory is propagated according to the (modified) FSSH algorithm with an ansatz for the rescaling direction as described above. For a particle moving in the 2-D plane, the magnetic forces are of the following form:

[TABLE]

For Method #1, the rescaling direction is

[TABLE]

For Method #2, we would ideally like to choose the direction $\bm{\nabla}\theta$, i.e. the $x$-direction, which we presume is the classical reaction coordinate. Unfortunately, with an arbitrary phase possible when delineating eigenstates, and without the knowledge of a global potential energy surface, isolating $\bm{\nabla}\theta$ is non-trivial. In the present case (for a general $\bm{d}$, see Appendix B), the vector norm $f(\eta)=||Re(e^{i\eta}\bm{d}_{01})||^{2}$ becomes

[TABLE]

Maximizing the above expression using $\bm{\nabla}\theta\cdot\bm{\nabla}\phi=0$, Method #2 chooses the rescaling direction to be:

[TABLE]

For most parameters below (except Fig. 7), we will usually operate in the regime whereby $(\partial_{x}\theta)^{2}>(\sin{\theta}\partial_{y}\phi)^{2}$, and so $\bm{n}_{01}$ will be in the $x$-direction. For our other parameters, we choose $B=3.0$, $\bm{r}_{init}=(-3,0)$, $\sigma=1.0$.

IV Results

We begin by investigating scattering processes where the average incoming momentum is along the $x$-direction: $\bm{p}_{init}=(p^{x}_{init},0)$. Because we initialize all dynamics to begin in the $x$-direction, we can learn about Berry’s phase effects by monitoring all dynamics in the $y$-direction.

In Fig. 2, we plot the final population and average momentum along the $x$-direction and the $y$-direction for each diabat as a function of initial $x$ momentum, $p^{x}_{init}$. As far as electronic populations are concerned, the case $A=0.02$ would appear to be in the diabatic regime at large velocities, where a significant percentage of trajectories stay on the initial diabat (diabat $1\rightarrow 1$); the case $A=0.1$ would appear to be in the adiabatic regime, where most trajectories stay on the initial adiabat (diabat $1\rightarrow 0$). The exact quantum dynamics results give a simple interpretation of Berry’s phase effect: motion on a given diabat (i.e. a switch of adiabat) does not lead to a finite momentum in the $y$-direction. By contrast, motion on two different diabats (no switch of adiabats) does lead to a finite momentum. All momentum changes are identical (with a value of $W=0.5$), irrespective of the values of $A$ and $p^{x}_{init}$.

Let us now turn to FSSH. For the population distribution and the $x$-direction momentum, we find that the modified surface hopping algorithm does capture the correct results (for both Methods #1 and #2). However, for the $y$-direction momentum, while the FSSH result is good for the case of small $A$, its error increases when $A$ is tuned to larger values. Again, there is no significant difference between the two rescaling ansätze.

Next, in Fig. 3, we investigate the same case but with a larger $W$: We set $W=5.0$ and plot the same observables as in Fig. 2. As far as the accuracy of FSSH is concerned, our conclusions are the same as for the case of $W=0.5$. Both FSSH ansätze capture the correct $p^{x}$ and the approximately correct $p^{y}$. The error increases as $A$ increases. For this model problem, the exact momentum change in the $y$-direction (for wave packets that do not switch adiabats) is again equal to $W$, only now $W=5.0$.

Finally, for a meaningful comparison of the two rescaling ansätze and as a means of differentiation, we turn to an alternative set of initial momentum conditions: All trajectories are initialized with a momentum $\bm{p}_{init}=(p^{x}_{init},-p^{x}_{init})$. As plotted in Fig. 4, the results of two ansätze become different: Method #2 almost captures the correct momentum distribution, while Method #1 consistently underestimates $p^{x}$ and overestimates $p^{y}$ on diabat 1. From this observation, we empirically infer that, if surface hopping is applicable with complex Hamiltonians, Method #2 must be more physically meaningful than Method #1. Evidently, the optimal rescaling direction is Method #2, which does not depend on any dynamical information.

V Discussion

V.1 Understanding Berry’s Magnetic Force: The case of small or moderately sized $W$

To incorporate Berry’s phase effects into semiclassical dynamics, we have used the well-known magnetic force ansatz Shankar (2012) in Eq. 4. To better understand this force in the context of semiclassical dynamics, note that, in Figs. 2 and 3, one finds that the exact momentum change in the $y$-direction is independent of $A$ and $p^{x}_{init}$. In fact, as long as $W$ is not too large, if the system is initialized on the upper surface, one always ends up with $p^{y}_{1}=W$ and $p^{y}_{0}=0$; if the system is initialized on the lower surface, one always ends up with $p^{y}_{1}=0$ and $p^{y}_{0}=-W$. Neither $A$ or $p^{x}_{init}$ has an effect on the $p^{y}$ results. These features are completely consistent with the fact that Berry’s phase is a topological (rather than dynamic) effect. Nevertheless, in Figs. 2 and 3, one finds that Berry’s topological phase has clear dynamic consequences.

Within the context of semiclassical dynamics, the ansatz of an extra magnetic force can partially handle these effects: Given the expression for $\bm{F}^{mag}_{1}$ in Eq. 14, if a trajectory is initialized on the upper surface and propagated adiabatically, i.e. without any hopping and assuming full transmission, we find that the final momentum in the $y$-direction is:

[TABLE]

Similarly, if a trajectory is initialized on the lower surface and propagated without hopping, we recover

[TABLE]

Clearly, Eqs. 18 and 19 are effectively the correct, semiclassical adiabatic limits; our ansatz for including Berry’s forces within FSSH appears reasonable.

Let us next address the question of whether semiclassical FSSH is trustworthy in practice in the limit of finite (or nonzero) hopping probabilities. In Fig. 2, in the case of a small $W$ value, we saw that FSSH almost captures the correct results but the agreement is not perfect. As we will show now, this non-agreement can be traced back to the very basic concept of independent FSSH trajectories with variable hopping positions. Consider for a moment the early surface hopping proposal by Tully and Preston Tully and Preston (1971), whereby a trajectory hops between adiabats only at a crossing point in the spirit of Laudau-Zener transition. In this case, it is quite easy to see that surface hopping should be nearly exact. On the one hand, for a trajectory that does not hop at the crossing point, the ending $p^{y}$ will be exactly $W$ given the limit of zero hopping (see Eq. 18). On the other hand, for a trajectory that hops at exactly the crossing point, half of the transmitted trajectory will run on one adiabat and half will run on the other adiabat. Thus, by symmetry of the Berry’s force, i.e., the fact that $\bm{F}^{mag}_{0}=-\bm{F}^{mag}_{1}$, the final $p^{y}$ will be [math]. Therefore, Tully-Preston surface hopping must be accurate for incorporating geometric phase, and any deviations in the FSSH the $p^{y}$ results must be caused by the fact that Tully’s FSSH algorithm allows trajectories to hop up and down, back and forth, multiple times in the coupling region; this complicated hopping picture no longer guarantees that the $y$-momentum induced by the Berry magnetic force will be accurate. In the end, the small inaccuracies in Figs. 2 and 3 appear inevitable if one sticks with the independent FSSH algorithm, even in the limit of small $W$.

V.2 The Limitations of the Modified FSSH

Next, let us consider larger $W$ values and/or non-perpendicular incoming velocities (so that $\bm{F}_{1}^{mag}$ is negative in the $x$-direction), where another feature can also appear: Reflection. Even though the adiabats are entirely flat, it is possible to observe reflection!

In Figs. 5 and 6, we let $p^{y}_{init}=p^{x}_{init}$ and we investigate both the transmitted and reflected particles, respectively. We find that both the exact quantum solution and the modified FSSH algorithms predict some amount of reflection provided that we apply the correct magnetic force in our FSSH algorithm. That being said, although Method #2 is still better than Method #1, neither method can fully capture the correct population and momentum quantitatively even when $A$ is small.

Finally, let us address the case of very large $W$.

In Fig. 7, we plot simulation results for $W=15$. For this case, an important nuance arises regarding to our FSSH algorithm. Unlike the case of small or medium $W$, where Method #2 is equivalent to rescaling in the $x$-direction, for the case of large $W$, Method #2 can actually rescale momenta along the $x$-direction for some coordinates but along the $y$-direction for others (see Eq. 17). As a means of assessing this unusual ansatz, we will introduce yet another rescaling scheme: Simple rescaling along the $x$-direction after a surface hop.

From the results in Fig. 7, we find that, when $W$ is large, no modified FSSH algorithm works well.333 Method #1 also fails (not shown). One is not even able to capture the electronic state populations as a function of $p^{x}_{init}$. One can conceive of two possible explanations for this dramatic failure: $(i)$ When $W$ is large, the complex Hamiltonian matrix oscillates rapidly with frequency $W$ as a function of the coordinate $y$, and so the dynamics may be outside the classical region, and quantum effects may be essential, as in the case of a time-dependent Hamiltonian with large frequency $\omega$. $(ii)$ It is also possible that we have not yet found the optimal approach for velocity rescaling after a hop. Understanding how and why FSSH fails in the case of large $W$ deserves further investigation.

V.3 Time Dynamics

Before concluding, let us turn to time dynamics rather than scattering probabilities. So far in this manuscript, we have focused on the asymptotic states after a scattering event – rather than the time dynamics of the underlying wave function during the scattering event.

To better understand the dynamics, in Fig. 8 we plot the time evolution for the populations and momenta on the diabats. The initial momentum $\bm{p}_{init}$ is set to be $(20,0)$. To obtain FSSH statistics on the diabatic surfaces, we use method 3 from Ref. Landry et al., 2013. This conversion method leads to intense oscillations if the wavefunction is initialized entirely on the upper diabatic surface, and to avoid such a numerical issue, we initialize the wavefunction with a slight superposition state: $99.5\%$ of the population is initialized on the upper diabatic surface, while $0.5\%$ are initialized on the lower one. From the data in Fig. 8, we find that, despite the fact that the scattering process is dynamically complicated, FSSH dynamics are not actually that bad (just as for real Hamiltonians): The population dynamics predicted by FSSH are reasonably accurate, and the overall trend of momentum dynamics are basically in agreement with the exact dynamics.

Lastly, let us turn to the adiabatic representation. To generate exact adiabatic momenta, we rotate the electronic wavefunctions from the diabatic representation to the adiabatic representation, using the analytical eigenvectors in Eq. 12. Note that the quantities $\left\langle\psi_{0}\middle|\bm{\nabla}\middle|\psi_{0}\right\rangle$ and $\left\langle\psi_{1}\middle|\bm{\nabla}\middle|\psi_{1}\right\rangle$ are usually non-zero (and of course purely imaginary). The contribution of these terms must be included when evaluating the momentum on the adiabatic surfaces. Thus, if the exact wavefunction is $\left|\Psi\right\rangle=C_{0}\left|\psi_{0}\right\rangle+C_{1}\left|\psi_{1}\right\rangle$, we estimate the exact momentum on the upper adiabatic surface to be

[TABLE]

Next, let us consider FSSH. Normally, because FSSH is defined in the adiabatic basis, one would expect FSSH to be most accurate in this representation. For FSSH, the adiabatic momenta are computed simply by averaging the momentum of all trajectories on a given adiabatic surface.

The dynamics on the adiabatic surfaces are plotted in Fig. 9. For the adiabatic populations, FSSH again captures accurate dynamics. For the momentum, however, the FSSH result on adiabat 0 is extremely inaccurate for early times. In theory, this error could arise because, at early times, the details of one wave packet spreading from one adiabat to another must reflect the quantum nature of matter waves. A simpler and more likely explanation, however, is that FSSH fails here simply because semiclassical dynamics treat $\bm{p}$ classically whereas exact quantum dynamics interprets momentum as a phase change (that can more naturally account for the presence of geometric phase). Either way, it is quite surprising that, physical observables (as calculated by FSSH) in a diabatic basis appear more accurate than those in an adiabatic basis.

VI Summary

To summarize, we have proposed a modified version of FSSH to incorporate non-adiabatic semiclassical systems with complex Hamiltonians. For a chemistry audience accustomed to non-adiabatic transitions, we have shown how to include complex Hamiltonians and Berry’s forces; for a physics audience accustomed to adiabatic dynamics with complex Hamiltonians, we have shown one means to take the non-adiabatic limit and including hopping. For motion along adiabatic surfaces, we invoke the usual concept of adiabatic “magnetic forces” to account for Berry’s phase (Eq. 6), and some evidence has been provided that this approach is compatible with standard FSSH. 444 Interestingly, we note that, even though FSSH is grounded in the notion of dynamics along adiabats, a few researchers have designed surface hopping schemes in a diabatic basis.Hack et al. (2000); Wang et al. (2015) Within such a diabatic framework, one would not be able to use “adiabatic magnetic forces” to account for Berry’s phase (as we have done here). Instead, one would need to account for Berry’s phase when adjusting velocities after a hop, and there is no guarantee that such an approach would be robust. ,Hack et al. (2000); Wang et al. (2015) For the momentum rescaling scheme, we compare three potential ansätze and show that Method #2 is the best rescaling scheme: after a hop $j\rightarrow k$, the momentum should be adjusted in the direction $Re(e^{i\eta}\bm{d}_{jk})$ where $\eta$ is chosen to maximize $||Re(e^{i\eta}\bm{d}_{jk})||^{2}$, which is the same direction as $\bm{\nabla}\theta$ for a two-state model. 555 Note that, if $W$ is large enough, as discussed in the context of Fig. 7, Method #2 and $\bm{\nabla}\theta$ may give different rescaling directions. Nevertheless, for such large $W$ values, FSSH does not appear to be accurate – again, see Fig. 7 – and so Method #2 would appear to be a robust ansatz that should be applicable for ab initio calculations.

Evidently, choosing a dynamical rescaling direction is not appropriate.

With these adjustments, our overall conclusion is that, a modified FSSH algorithm can capture many important non-adiabatic dynamical features (e.g. the scattering probabilities and the approximate scattering momenta), but FSSH cannot capture a few features (e.g. the detailed early time dynamics of momentum transfer).

VII Open Questions

With the above summary in mind, several questions now present themselves. On the practical side, the first methodological question one must pose is: Have we constructed the optimal FSSH algorithm or is there another, better option available for the case of a complex electronic Hamiltonian? Considering the errors in the $y$-momentum in Figs. 2 and 3 and the discussion of independent trajectories in Sec. V.1, we note that Truhlar et al have constructed an FSSH algorithm with time uncertaintyJasper et al. (2002) which was designed to introduce a small amount of time non-locality. Would a similar approach help improve FSSH in this case and reduce the number of hops in the coupling region? Or is it simply impossible to model Berry’s phase well with independent trajectories, given that Berry’s phase is geometric and topological (and therefore intrinsically non-local)?

Second, again on the practical side, a modern FSSH implementation can avoid calculating derivative couplings unless a hop is required;Jain et al. (2016) as far as propagating time dependent Schödinger equation, $\bm{d}\cdot\bm{p}/m$ is enough. Unfortunately, in the case of a complex Hamiltonian with Berry’s forces, apparently one must calculate $\bm{d}$ at every time step in order to evaluate $\bm{F}^{mag}_{j}$. One must wonder: Is there a practical and efficient approach to construct such a Berry force easily, ideally a scheme that will be stable with a large number of electronic states and will avoid the trivial crossing problem?Fernandez-Alberti et al. (2012); Nelson et al. (2013); Plasser et al. (2012); Wang and Prezhdo (2014); Meek and Levine (2014); Jain et al. (2016)

Third, on the theory side, one must also wonder: Can any of our proposed extensions of FSSH be tied back to a more rigorous theory of quantum mechanics? For the case of a real electronic Hamiltonian, our research group and the Kapral research group have successfully tied FSSH back to the QCLEKapral (2016); Subotnik et al. (2013). However, the QCLE is a first order expansion that cuts off at zeroth order in $\hbar$, whereas Berry’s phase requires a second-order expansion: Note that the magnetic force in Eq. 6 is first order in $\hbar$. Can we relate an extended version of FSSH to an extended version of the QCLE for the case of complex Hamiltonians?

Fourth, according to Figs. 5 and 6, the magnetic force in Eq. 6 can lead to wave packet separating as trajectories on different adiabatic surfaces are turned in different directions, some transmitted and some reflected. Thus, the sharp reader will no doubt isolate yet another question. Recall that, when deriving FSSH from the QCLE, the question of decoherence and wave packet separation arises naturally.Kapral (2016) After all, wave packets on different adiabatic surfaces feel different static, adiabatic forces that lead to separation eventually; and for years, many researchers have constructed practical solutions for incorporating decoherence into FSSH to account for such effects.Bittner and Rossky (1995); Schwartz et al. (1996); Fang and Hammes-Schiffer (1999b); Subotnik and Shenvi (2011); Wang et al. (2016); Hack and Truhlar (2001); Volobuev et al. (2000); Jasper and Truhlar (2005); Prezhdo and Rossky (1997b) For the present paper, however, we now see a new phenomenon: With Berry’s forces, wave packet separation is caused by wave packets on different surfaces feeling different magnetic forces that depend on velocity. Furthermore, these “magnetic” forces appear only in the strong coupling region, which negates our usual understanding of decoherence being a phenomenon that emerges after wave packets pass through coupling region and only thereafter move apart in different directions. Zhu et al. (2004) Thus, another immediate question is how should we appropriately model such magnetically induced decoherence within FSSH so as to recover the correct dynamics.

Given the inherent difficulties of including decoherence within FSSH, the questions above lead to a fifth question: Is it possible that a different mixed quantum classical scheme might strongly outperform FSSH for the case of complex Hamiltonians? In particular, for problems of decoherence, ab initio multiple spawning (AIMS) is a more natural ansatz.Martínez et al. (1996); Ben-Nun et al. (2000) And yet, AIMS is most efficient in an adiabatic basis, where single valued wave functions can be difficult to find. Interestingly, there has been a great deal of work investigating conical intersection’s geometric phase and choice of basis within AIMS for real Hamiltonians, and the overall conclusion appear to be that we should run dynamics with electronic wavefunctions chosen at a single location.Meek and Levine (2016) Thus, one can ask, can the results in Ref. Meek and Levine, 2016 for adiabatic AIMS be easily extended to work with complex Hamiltonians?

The final, sixth question is perhaps most exciting of all. On the experimental front, one must wonder: Can any of the dynamics predicted in Sec. IV above be detected experimentally? For instance, the numerical model above suggests that, whenever an electronic transition (in the $x$-direction) occurs between two electronic states with spin, one ought to find a signature of nuclear or vibrational motion (in the $y$-direction) as arising from Berry’s phase for the case of a molecule in a magnetic field – provided that the transition occurs in the normal regime where an electron changes character along a single adiabat. Vice versa, no such signature should be observable for a transition in the inverted regime where an electron changes character but the adiabat also changes. Can this dichotomy be seen experimentally? Can we find realistic molecular systems with large enough susceptibilities such that, in very large magnetic fields, we will observe dynamical Berry phase effects? Or, if we recall that Marcus theory assumes a threshold amount of nuclear friction, a pessimist must ask: Will the inevitable presence of some nuclear friction eliminate all such effects? And lastly, how will these features behave when the complex phase is more complicated, so that $\partial_{y}\phi$ is not a constant (as assumed above)? These fascinating experimental and theoretical questions will hopefully be answered in the near future.

Appendix A Method #3 Results

Here we briefly present scattering results for Method #3 in II.2.2.

In Fig. 10, we report results for the case $W=0.5$ with Methods #3 and #2. Clearly, the results indicate that neither the correct population nor the correct momentum can be captured by Method #3. Now, at first glance, it might seem that the $y$-momentum on diabat 0 is correctly captured. This instinct is merely an illusion, however, as this “accurate” results is caused only surreptitiously from the fact that, at the end of the simulation, the particles remaining on the upper adiabatic surface (i.e. diabat 0) are mostly those trajectories that never hop, and so the average $y$-momentum will always go to the correct answer (as induced by the magnetic force). This correct answer is the zero hopping limit (or adiabatic limit) as discussed in Sec. V.1 that arises from simple classical mechanics (and ignoring all surface hops).

Overall, even though it might appear natural, Method #3 does not coincide with the correct physical picture.

Appendix B Method #2 with a General Derivative Coupling

Here, we analyze Method #2 for a general, two-state diabatic problem. We denote a general (complex) derivative coupling vector as $\bm{d}\equiv\bm{d}_{R}+i\bm{d}_{I}$. When using Method #2, we maximize the following within the interval $\eta\in[0,\pi)$:

[TABLE]

Setting $f^{\prime}(\eta)=0$ tells us $\eta$ should satisfy

[TABLE]

Within the $[0,\pi)$ interval, there exist two solutions: $\eta_{0}$ and $\eta_{1}=\eta_{0}+\pi/2$, with $\eta_{0}\in[0,\pi/2)$. Using the second derivative $f^{\prime\prime}(\eta)<0$, we find that maximizing $f(\eta)$ requires that $\eta$ must satisfy

[TABLE]

Thus, for a general $\bm{d}$, the solution for $\eta$ should satisfy both Eq. 22 and Eq. 23.