Restoring electronic coherence/decoherence for a trajectory-based nonadiabatic molecular dynamics

By utilizing the time-independent semiclassical phase integral, we obtained modified coupled time-dependent Schrödinger equations that restore coherences and induce decoherences within original simple trajectory-based nonadiabatic molecular dynamic algorithms. Nonadiabatic transition probabilities simulated from both Tully’s fewest switches and semiclassical Ehrenfest algorithms follow exact quantum electronic oscillations and amplitudes for three out of the four well-known model systems. Within the present theory, nonadiabatic transitions estimated from statistical ensemble of trajectories accurately follow those of the modified electronic wave functions. The present theory can be immediately applied to the molecular dynamic simulations of photochemical and photophysical processes involving electronic excited states.

where R is an N-dimensional vector of nuclear coordinates and an overdot denotes a time derivative (N-dimensional velocity) with the sum over all adiabatic electronic states. Tully's fewest switches 1 and the semiclassical Ehrenfest 2 algorithms, which are the two representative methods, provide a powerful tool to perform nonadiabatic molecular dynamics simulations with simple trajectory-based approaches. These trajectory-based nonadiabatic molecular dynamics methods with various modified versions have been successfully applied to photochemical and photophysical related molecular spectra and reaction dynamics for large systems [3][4][5][6][7][8][9][10][11][12][13][14][15][16][17] . Tully's fewest switches method propagates nonadiabatic trajectory on each adiabatic electronic state with trajectory surface hopping from one adiabatic potential energy surface to another. The semiclassical Ehrenfest method propagates nonadiabatic trajectory on averaged adiabatic electronic states in which a single mean-field potential energy surface governs nuclear motion. Both methods suffer an overcoherence problem in electronic wavefunction (or density matrix ρ kj in eq. (1)) when nonadiabatic trajectory passes through non-zero region of nonadiabatic coupling vector. In order to reproduce exact coherent motion of quantum wavefunction propagating on multiple adiabatic electronic states, the semiclassical density matrix described by eq. (1) should decohere. Numerous algorithmic coherence/decoherence schemes in the literature have been designed for nonadiabatic trajectory coupled with electronic motion [18][19][20][21][22][23][24][25][26][27][28][29][30][31][32] . All of these approaches can improve decoherence/coherence effects of the electronic wavefunction to some extent but at either high computational cost or through the use of complicated algorithms. in a unified way for both Tully's fewest switches and the semiclassical Ehrenfest methods. Let us exam solution of eq. (1) in the regions where nonadiabatic coupling vector is equal to zero, we have kj kj t k j 0  Let us compare the phase integral in eq. (2) with the conventional time-independent Jeffreys-Wentzel-Kramers-Brillouin 33 where μ is reduced mass, E is total energy and T = (E − U) is kinetic energy. From trajectory-based approach, R in eq. (3) can be considered as a classical trajectory propagating along a one-dimensional curved space. Alternatively, we can obtain the same relation as in which U(R) is an effective potential energy surface for nonadiabatic trajectory. For Tully's fewest switches approach, U(R) is a single adiabatic potential energy surface U j (R) or U k (R) on which trajectory is propagating, and for semiclassical Ehrenfest approach, U(R) is defined as an average potential energy surface Finally, we obtain modified coupled time-dependent Schrödinger equations,  (4) shows that time-dependent and time-independent approaches have the same limit in the high-energy regime. We can expect that both original and modified coupled Schrödinger equations should agree each other for the high-energy regime, namely E − U ≫ |U k (R) − U j (R)|).
Using a similar approach to the semicalssical Ehrenfest method, the symmetrical windowing quasi-classical approach 27 (that is as simple as the present theory) restores coherence/decoherence pretty well with windowing parameter γ = 0.366. However, numerical tests in Supplementary Note 1 show that the present modified Ehrenfest approach is slightly better than the symmetrical windowing quasi-classical approach. Using a similar approach to the Tully's fewest switches trajectory surface hopping method, the Gaussian wavepackets phase correlation method 28 (that is more complicated than the present theory) restores coherence/decoherence in a similar way as the present one, and actually two are the same in one-dimensional case. However, the two are different in multidimensional case, their diagonal element in Hamiltonian does not approach limit U k (R)−U j (R) in the high-energy regime, and their decoherent term (p 1 · p 2 /m) requires calculation of classical momentum on two adiabatic potential energy surfaces simultaneously. Numerical tests in Supplementary Note 3 show that the present modified the Tully's fewest switches approach can nicely reproduce quantum oscillation for the certain two-dimensional model system.
We select four well-known model problems in the following to perform both semiclassical Ehrenfest and Tully's fewest switches calculations from original eq. (1) and modified eq. (6). Model 1 describes the electronic transitions in the non-crossing case of adiabatic potential energy surfaces initially developed by Rosen and Zener 34 . Model 2 and Model 3 describe the electronic transitions in the avoided-crossing case of adiabatic potential energy surfaces initially developed by Landau 35 , Zener 36 and Stückelberg 37 . Model 4 describes the electronic transitions in the crossing case with the peculiar degeneracy of adiabatic potential energy surfaces initially developed by Renner 37,38 . All these pioneer studies focused on developing analytical formula for nonadiabatic transition probability from various mathematical methods nicely documented by Child 39 . The present study focuses on simulating nonadiabatic transition probability for all cases numerically solving original eq. (1) and modified eq. (6).
We compute the overall nonadiabatic transition probability defined as the probability of starting on the lower adiabatic potential energy surface at x = −∞ and finishing on the upper adiabatic potential energy surface at x = + ∞ . Accurate quantum mechanical calculations for the four one-dimensional two-state models are performed using the conventional time independent close-coupling method. Semiclassical calculations are performed by using the fourth-order Runge-Kutta method for numerically integrating the trajectories as well as coupled time-dependent Schrödinger equations. For a given total energy, the Tully's fewest switches method  with the parameters chosen to be A = B = 0.025 Hartree, x 0 = 3.0 Bohr, and C = 0.7 Bohr −2 . This model can induce electronic coherence from overlapping of two semiclassical wavepackets in between two non-crossing peaks at x = ± 3.0 Bohr as shown in Fig. 1(a). The results simulated by Tully's fewest switches (see Fig. 2(a)) and semiclassical Ehrenfest (see Fig. 2(b)) within the modified semiclassical eq. (6) follow exact quantum oscillations and amplitudes of the overall nonadiabatic transition probabilities very well. The results simulated from the original semiclassical eq. (1) cannot reproduce exact quantum results, and besides Tully's fewest switches and semiclassical Ehrenfest methods do not agree each other for oscillations and amplitudes of the overall nonadiabatic transition probabilities.

Dual Landau-Zener-Stückelberg avoided-crossings. Model 2 is defined by two diabatic potentials
having dual crossings given in the diabatic representation: This model can induce different electronic coherence from overlapping of two semiclassical wavepackets in between two avoided-crossings at x = ± 2.07 Bohr as shown in Fig. 1(b). The results simulated by Tully's fewest switches (see Fig. 3(a)) and semiclassical Ehrenfest (see Fig. 3(b)) within the modified semiclassical eq. (6) follow exact quantum oscillations and amplitudes of the overall nonadiabatic transition probabilities very well. The results simulated from the original semiclassical eq. (1) cannot reproduce exact quantum results, and besides Tully's fewest switches and semiclassical Ehrenfest methods show very different oscillations and amplitudes of the overall nonadiabatic transition probabilities. This model system has been well studied in the   Simple avoided crossing. Model 3 is defined by two diabatic potentials having one simple crossing given in the diabatic representation: with A = C = 0.01 Hartree, B = 1.6 Bohr −1 , and D = 1.0 Bohr −2 . One simple diabatic crossing is appeared at x = 0 as shown in Fig. 1(c). However, the strong Gaussian-type diabatic coupling can produce two vague avoided crossings around x = ± 1.0 Bohr, and thus it can induce small electronic coherence at low energies as shown in Fig. 4 The results simulated by Tully's fewest switches (see Fig. 4(a)) and semiclassical Ehrenfest (see Fig. 4(b)) within the modified semiclassical eq. (6) can well reproduce this small oscillation and amplitudes of nonadiabatic transition probabilities. The results simulated from the original semiclassical eq. (1) cannot well reproduce this small oscillation, although Tully's fewest switches and semiclassical Ehrenfest methods almost agree with the exact quantum amplitudes of the overall nonadiabatic transition probabilities for high energies, but two are not exactly the same.
Renner-Teller crossing. Model 4 is defined by two diabatic potentials having same as V 11 (x) and V 22 (x) in eq. (9). However, diabatic coupling is changed as  modified methods cannot follow exact calculation of the overall nonadiabatic transition probability as shown in Fig. 5. This means that the restoring term in eq. (6) is still not good enough for describing electronic transitions in significant degenerate case and this term basically comes from approximation of kinetic energy operator in nuclear degree of freedom. Further study should be carried out in the near future.
Concluding remarks. In three out of the four model systems given above, we have shown that the results simulated from the modified semiclassical eq. (6) follow exact electronic coherence as well as amplitude of the overall nonadiabatic transition probabilities, and both the semiclassical Ehrenfest and Tully's fewest switches methods agree with each other and even for small oscillation in Model 3. On the other hand, the results simulated from the original semiclassical eq. (1) cannot follow exact electronic coherence, and besides the semiclassical Ehrenfest (Tully's fewest switches) method shows slight greater (smaller) amplitude of the overall nonadiabatic transition probabilities than exact quantum results. This can be seen clearly from Model 3 in the region of monotonically increasing of the overall nonadiabatic transition probabilities against energy. More tests have been performed with changes of potential parameters, and conclusion is the same for the modified coupled Schrödinger equations. The present modified eq. (6) only modifies diagonal element by in eq. (1) while preserving original simplicity of Tully's fewest switches and semiclassical Ehrenfest algorithms. For instance, detailed balance behavior in the present Tully's fewest switches should follow Tully's fewest switches and global flux surface hopping 41 , while the present semiclassical Ehrenfest should follow detailed balance of symmetrical quasi-classical treatment of the Meyer-Miller (MM) model 42 . In the multi-state nonadiabatic molecular dynamic simulation, the global flux surface hopping 43 shows promising accuracy in comparing Tully's state-to-state surface hopping. This can be immediately applied to the present coupled electronic Schrödinger eq. (6) to perform multi-state nonadiabatic molecular dynamic simulation. However, it should be noticed that when two adiabatic potential energy surfaces significantly degenerate at crossing zone with Renner-Teller coupling, the present modified method still fails. Otherwise, we conclude that in the present theory both semiclassical Ehrenfest and Tully's fewest switches algorithms are shown to work equally well and the both follow electronic coherence/decoherence the same as exact quantum results. The present theory can be immediately applied to nonadiabatic molecular dynamic simulation for photochemical and photophysical processes involving with electronic excited states.