We propose an approximate method for evaluating the importance of non-Born–Oppenheimer effects on the quantum dynamics of nuclei. The method uses a generalization of the dephasing representation (DR) of quantum fidelity to several diabatic potential energy surfaces and its computational cost is the cost of dynamics of a classical phase space distribution. It can be implemented easily into any molecular dynamics program and also can utilize on-the-fly *ab initio* electronic structure information. We test the methodology on three model problems introduced by Tully and on the photodissociation of NaI. The results show that for dynamics close to the diabatic limit, the decay of fidelity due to nondiabatic effects is described accurately by the DR. In this regime, unlike the mixed quantum-classical methods such as surface hopping or Ehrenfest dynamics, the DR can capture more subtle quantum effects than the population transfer between potential energy surfaces. Hence we propose using the DR to estimate the *dynamical* importance of diabatic, spin-orbit, or other couplings between potential energy surfaces. The acquired information can help reduce the complexity of a studied system without affecting the accuracy of the quantum simulation.

*Introduction.* The nonadiabatic effects play an important role in many chemical phenomena and often must be taken into account in accurate calculations of molecular properties such as spectra or reaction rates.^{1} Many nonadiabatic quantum (QM) simulations are performed in the diabatic (or quasidiabatic) basis, which offers several computational advantages, especially in the vicinity of conical intersections of adiabatic potential energy surfaces (PESs). Usually, more PESs must be included to describe a system of interest accurately. One may ask: Which diabatic surfaces are important? How much is the result affected by neglecting the less important surfaces? The method we propose below quantifies the importance of couplings between PESs and thus can help in answering these questions.

A direct way to determine whether the coupling of diabatic PESs affects the result of a simulation is to compute the desired quantity by running QM dynamics for both the uncoupled and coupled systems and to compare the results. Here, we propose a more general approach, which can give the information about the importance of couplings for all observables at the same time. It consists in comparing the wave functions by computing the QM fidelity, defined as^{2}

In the general setting, $\psi 0(t)$ and $\psi p(t)$ are wave functions evolved in the unperturbed and perturbed systems, respectively. In our case, the “unperturbed” represents the *uncoupled* system, “perturbed” the *coupled* system, and $\psi $ denotes the full molecular wave function, i.e., it includes both nuclear and electronic degrees of freedom. Values of $FQM(t)$ close to unity imply that the perturbation is not important: The uncoupled Hamiltonian can be used in quantitative simulations. Values significantly below unity imply that the perturbation is important: The affected PESs and couplings should be included in the simulation.

Instead of computing fidelity directly from the definition (1), one can use the dephasing representation (DR),^{3–5} which is an efficient semiclassical approximation of fidelity.^{6,7} Recently, the DR was used successfully to evaluate the accuracy of QM dynamics on a single but approximate PES.^{8,9} Below, we generalize this method to several surfaces and use it to evaluate how the dynamics is affected by couplings between the diabatic PESs.

Unlike the computational cost of a direct QM calculation, the cost of the DR does not grow exponentially with the number $D$ of degrees of freedom.^{5,8} Hence, the DR can be used for many-dimensional systems inaccessible to current methods of QM dynamics. An advantage of the DR in comparison with mixed QM-classical methods for nonadiabatic dynamics, such as the Ehrenfest dynamics, various surface hopping methods,^{10,11} or methods in which the classical limit is obtained by the linearization of the QM propagator,^{12} is that the DR, being a semiclassical method, takes the nuclear coherence effects into account at least approximately. Other semiclassical approaches to nonadiabatic dynamics exist, including, among others, the multiple-spawning methods^{13} or methods^{14} based on the Miller–Meyer–Stock–Thoss classical electron model. The advantage of the DR is that unlike the majority of semiclassical approaches, the DR does not require the Hessian of the potential energy, which is the most expensive element of first-principles semiclassical dynamics methods (see, e.g., Ref. 15).

Clearly, the above-mentioned advantages do not come for free. First, unlike most other approaches noted above, the DR is not a general dynamical method. It can only describe properties expressible in terms of fidelity. Second, the DR in the diabatic basis is expected to work accurately mainly when the coupling is relatively weak or the nuclear velocities are large, i.e., when the dynamics is close to the diabatic limit. Nevertheless, in many cases, the DR works very well even for strong perturbations that lead to completely different classical phase space distributions for uncoupled and coupled Hamiltonians.

Applications for which the DR method is particularly well suited include the chemical reactions which proceed near the diabatic limit. The photodissociation of bromopropionyl chloride or bromoacetyl chloride^{16} serve as examples. Another possible application would be to quantify the relative importance of Hamiltonian coupling terms of various origins, e.g., of the spin-orbit versus diabatic couplings. For instance, in the $Cl(P2)+H2$ reaction, the spin-orbit coupling terms clearly dominate over diabatic couplings, which can therefore be neglected in a simulation.^{17}

*Theory.* The starting point is the Hamiltonian describing nuclear motion in a molecule, expressed in the diabatic basis. To compute the decay of fidelity due to the coupling between $n$ PESs, the Hamiltonian $n\xd7n$ matrix is split into the uncoupled (i.e., diagonal) part $H\u0302diag$ and the coupling (i.e., off-diagonal) part $\Delta V\u0302$,

In Eq. (3), $T\u0302$ is the diagonal nuclear kinetic energy matrix and $V\u0302diag$ contains the diabatic PESs, which are uncoupled in $H\u0302diag$ and coupled in $H\u0302$ by the elements of $V\u0302offdiag$. (Bold face denotes $n\xd7n$ matrices; hat $\u2009\u0302$ denotes operators.)

To derive the DR, one starts from the expression for QM fidelity amplitude, applicable to both pure and mixed states^{5} and generalized to the multi-PES setting

where $\rho \u0302$ is the density operator of the initial state. Generalizing the derivation from Ref. 18, fidelity amplitude can be written exactly as

where $\rho W$ is the Wigner transform of the initial state $\rho \u0302$

and $x$ denotes the point $(q,p)$ in the $2\xd7D$-dimensional phase space. Approximating the Wigner transform of the product of the time-evolution operators, one arrives at the DR expression

where $\Delta S(x0,t)$ is the action at time $t$ due to Wigner representation $\Delta VW$ of $\Delta V\u0302$ along the trajectory $x\tau $ of $Hdiag$. Note that $\Delta VW$, $\Delta S$, and $\rho W$ are matrix quantities. If $\Delta V\u0302$ contains only diabatic coupling elements, then $\Delta V\u0302\u2261\Delta V(q\u0302)$ and $\Delta VW(x)=\Delta V(q)$. For initial Gaussian wave packets, the Wigner function equals the classical phase space density which is strictly a probability distribution.

Since $Tr\u222bdx0\rho W(x0)=1$, it follows from Eq. (6) that $fDR$ can be computed as a Monte Carlo average $\u27e8exp(\u2212i\Delta S(x0,t)/\u210f)\u27e9\rho W(x0)$, with initial conditions sampled from the Wigner distribution $\rho W(x0)$ of the initial state. However, as $\Delta S$ in Eq. (7) is much smaller than the action in other semiclassical methods, the DR alleviates the notorious “sign problem.”

If the initial state lies on a single surface $n$, i.e., if $[\rho W]nj=0$ and $[\rho W]jj=0$ for $j\u2260n$, then no other elements of $V=Vdiag+Voffdiag$ than $Vnn$ and $Vnj,j\u2260n$ enter the calculation. This means that no information about other diagonal elements $Vjj$ is used in the DR calculation. Thus, except for special cases, $FDR$ is expected to approximate $FQM$ accurately only when the detailed structure of the remaining PESs does not significantly affect the dynamics on $Vnn$.

*Results and Discussion: Tully's model problems.* First, the method is tested on three one-dimensional model potentials proposed by Tully^{11} to cover the most important characteristics of nonadiabatic transitions. Diabatic and adiabatic PESs as well as the coupling terms $V12$ for Tully’s problems A and C are shown in Figs. 1 and 2; further details can be found in Ref. 11.

In all cases, the initial wave packet was a Gaussian with approximately 10% dispersion in momentum and located on the lower (problems A and B) or upper (problem C) energy diabatic PES in the asymptotic region without diabatic couplings. The equations of motion were integrated until the wave packet left the interaction region.

The simple avoided crossing model (problem A in the original paper^{11}) represents the most often encountered situation. As can be seen in Fig. 1, the survival probability $P1,QM=Tr\rho \u030211(t)$ on the initial PES is nearly equal to the QM fidelity, suggesting that the fidelity decay is caused almost exclusively by the transition to the second diabatic PES. Close to the diabatic limit, $FDR$ accurately reproduces $FQM$. In practice, the error of $FDR$ comes from the intrinsic error of the approximation and from the statistical error due to the finite number of trajectories. In this case, the intrinsic relative error is $\u223c0.5%$ when fidelity decays by 10% to 0.9, and even less for the wave packets with higher initial momentum and hence fidelity. The relative statistical error remains below 1% even when only a single trajectory is used, since most of the fidelity decay is due to the transitions to the second surface. For slower wave packets, fidelity decreases, and the agreement between $FQM$ and $FDR$ stays qualitative, with $FDR$ always decaying faster than $FQM$. Finally, for a fixed initial momentum, decreasing the coupling $V12$ results in a slower decay of fidelity but the relative error of $1\u2212FDR$, i.e. $(FQM\u2212FDR)/(1\u2212FQM)$ is approximately independent of the coupling. Conversely, if $V12$ is increased substantially, fidelity initially decays quickly to zero, which is captured well by the DR. Subsequently, fidelity may rise again, which is usually not well reflected by the DR.

In the dual avoided crossing model (problem B, not shown), the DR works best again for high energy wave packets. For those, the detailed structure of the second PES does not significantly affect the motion on the first PES. At low energies, where a significant transfer of probability density to the diabatic surface $V22$ and back occurs, the DR fails to reproduce the QM fidelity even qualitatively. Nevertheless, in this case, both $FDR$ and $FQM$ initially decay almost to zero (although they might rise again later), correctly reflecting that the coupling is important and should not be neglected.

A very interesting situation occurs in the extended coupling model (problem C), where the two diabatic PESs $V11$ and $V22$ are almost equal (they differ just by a small constant shift) but the adiabatic surfaces are well separated due to the coupling. At very low energies of the wave packet [Fig. 2(b)], fidelity $F$ goes to zero despite that the survival probability $P1$ on the upper surface remains close to unity. Therefore, in contrast to previously discussed cases, it is not possible to estimate the extent of the nondiabaticity of QM dynamics on the basis of electronic transitions only. Unlike survival probability, fidelity describes the nondiabaticity effects on the coherent nuclear dynamics correctly. It is remarkable that the DR describes $F$ so accurately despite that the DR dynamics on the diabatic surface ignores the reflection of a large part of the QM wave packet from the upper adiabatic surface. The decrease of fidelity is correctly captured by $FDR$ for all energies of the wave packet. At higher energies [Fig. 2(c)], fidelity converges to zero after several oscillations, whereas the survival probability converges to 1/2. At very high energies (not shown), two wave packets moving on the two PESs interfere for a long time; $FDR$ and $FQM$ agree and oscillate between zero and unity. Surface hopping and Ehrenfest dynamics might predict $P1$ quite well, but since $V11\u2248V22$, surface hopping would incorrectly predict that $F\u2248P1$. While Ehrenfest dynamics would predict a decay of $F$ below $P$, it is unlikely that it would capture $F$ quantitatively since the QM wave packet splits into faster and slower components, whereas Ehrenfest dynamics uses a single mean-field surface. One of the reasons why the DR performs so well in the extended coupling model is the similarity of the diabatic PESs in the coupling region. If the surfaces are different, e.g., when $V22=cx$, the fidelity decay at low energies is still well approximated. However, at higher energies, $FDR$ oscillations slowly dephase from $FQM$, since the DR neglects effects of $V22$.

*Photodissociation of NaI.* We also applied the methodology to the photodissociation of NaI using a two-surface model of Engel and Metiu^{19} (see Fig. 3). In the original experiment of Mokhtari *et al.*,^{20} the molecule was excited by the light in the 310–390 nm range, which led to only a weakly nonadiabatic motion of the wave packet on the excited surface. In this regime (corresponding to$T0\u22480.07\u2002a.u.$ at the plateau of the PES), far from the diabatic limit, $FDR$ does not quantitatively reproduce $FQM$. However, when the momentum is approximately ten times higher and the dissociation of NaI is almost diabatic, $FDR$ reproduces $FQM$ rather well [Fig. 3(b) and 3(c)].

*Computational details.* QM dynamics calculations used the first-order split operator method. Classical trajectories were computed using the first-order symplectic Euler algorithm. The time step varied from 0.02 to 1 a.u., depending on the system and initial momentum.

To guarantee convergence, 16 384 classical paths were used to produce the DR plots. The statistical error of $fDR$ for a two-surface system due to a finite number $N$ of paths is given by

where $\phi =\Delta S12(x0,t)/\u210f$ is the phase accumulated along a trajectory and $\u27e8cos(\phi )\u27e9=fDR$. As demonstrated on Tully’s model A, where a single trajectory was sufficient, a much lower number of trajectories than 16 384 is needed to obtain an accurate estimate of fidelity in cases where fidelity stays close to unity. Equation (8) implies that for a given value of $\sigma stat$ and $fDR$, the number $N$ of trajectories needed is approximately independent of dimensionality.

*Conclusions.* Presented results demonstrate the utility of the DR in analyzing the molecular QM dynamics involving multiple PESs. On one hand, in the nearly diabatic regime, $FDR$ accurately approximates $FQM$. On the other hand, in systems far from the diabatic limit, $FDR$ decays quickly to zero and thus detects the importance of nondiabatic couplings although it may not reproduce $FQM$ accurately. Hence, the method can be used to establish the level of nondiabaticity of QM dynamics, without the need for a QM dynamics simulation. In fact, we propose the condition $F\u22481$ (instead of the standard requirement of high survival probability) as the rigorous definition of the diabatic limit (see Fig. 2). Nevertheless, for single avoided crossings, e.g., fidelity could be used to estimate the survival probability and hence the branching ratios (see Figs. 1 and 3).

The DR calculation can be performed easily for all systems accessible to classical molecular dynamics and for which coupling elements $Vij,i\u2260j$ are available. However, it remains to be verified how the method will perform in higher-dimensional systems and in systems with more than two important surfaces. The first issue was partially addressed in Ref. 8 where the DR was applied to the two-dimensional photodissociation of $CO2$, confirming that neither the accuracy nor the number of classical trajectories needed are significantly affected by increased dimensionality.

At the moment, we are exploring a related and complementary problem to that of the present paper, namely, whether the DR in the adiabatic basis can estimate the level of nonadiabaticity of QM dynamics near the adiabatic limit.

*Acknowledgments.* This research was supported by the Swiss NSF (Grant No. 200021_124936/1) and by the EPFL. We thank Ivano Tavernelli for helpful discussions.