Which form of the molecular Hamiltonian is the most suitable for simulating the nonadiabatic quantum dynamics at a conical intersection?

Choosing an appropriate representation of the molecular Hamiltonian is one of the challenges faced by simulations of the nonadiabatic quantum dynamics around a conical intersection. The adiabatic, exact quasidiabatic, and strictly diabatic representations are exact and unitary transforms of each other, whereas the approximate quasidiabatic Hamiltonian ignores the residual nonadiabatic couplings in the exact quasidiabatic Hamiltonian. A rigorous numerical comparison of the four different representations is difficult because of the exceptional nature of systems where the four representations can be defined exactly and the necessity of an exceedingly accurate numerical algorithm that avoids mixing numerical errors with errors due to the different forms of the Hamiltonian. Using the quadratic Jahn-Teller model and high-order geometric integrators, we are able to perform this comparison and find that only the rarely employed exact quasidiabatic Hamiltonian yields nearly identical results to the benchmark results of the strictly diabatic Hamiltonian, which is not available in general. In this Jahn-Teller model and with the same Fourier grid, the commonly employed approximate quasidiabatic Hamiltonian led to inaccurate wavepacket dynamics, while the Hamiltonian in the adiabatic basis was the least accurate, due to the singular nonadiabatic couplings at the conical intersection.

Unlike the adiabatic electronic states, which are only coupled through the nonadiabatic couplings, the quasidiabatic states have both the residual, presumably small, nonadiabatic couplings and the diabatic couplings-the off-diagonal elements of the potential energy matrix.However, the exact quasidiabatic Hamiltonian is rarely used.Instead, the approximate quasidiabatic Hamiltonian, which ignores the residual nonadiabatic couplings, is almost always used due to its simplicity.The separable form of the approximate quasidiabatic Hamiltonian allows using a wider range of time propagation schemes, [33][34][35][36] including the well-known split-operator algorithm, [37][38][39] but ignoring the residual nonadiabatic couplings decreases the accuracy. 40The strictly diabatic Hamiltonian, with only diabatic couplings and no nonadiabatic couplings, would be the most suitable for the nonadiabatic quantum dynamics simulation.However, in typical systems, the strictly diabatic states only exist in general when an infinite number of electronic states are considered. 41,42vantages and disadvantages of various Hamiltonians have been explored by numerous comparisons of the nonadiabatic quantum dynamics simulated with different Hamiltonians: To name a few, there exist comparisons between the strictly diabatic and approximate quasidiabatic Hamiltonians, [42][43][44][45][46] between the adiabatic and approximate quasidiabatic electronic states |n(Q) , obtained by solving the time-independent Schrödinger equation one can establish an approximate ansatz for the solution of the time-dependent Schrödinger equation here, ψ n (Q, t), V n (Q), and e iAn(Q) are the time-dependent nuclear wavefunction, potential energy surface, and coordinate-dependent phase factor associated with the nth adiabatic electronic state.The Born-Huang expansion 62 in Eq. ( 2) is not exact unless the sum includes an infinite number of terms but can be very accurate if a finite number S of electronic states are chosen wisely. 10,63,64e is free to choose an overall phase A n (Q) in Eq. ( 2) because if |n(Q) is a normalized solution of Eq. ( 1), then so is e iAn(Q) |n(Q) .2][23][24][25][26][27][28][29][30][31][32] In Sec.S4 of the supplementary material, we show that neglecting this double-valuedness of the wavepackets is detrimental to accuracy.Instead, in what follows, we set phases A n (Q) appropriately (see Sec. S1 of the supplementary material) to ensure the single-valuedness of both the adiabatic states and the wavepackets.In other words, we include the "geometric phase" in the adiabatic states in order to obtain the best possible results in the adiabatic representation.From now on, we absorb the overall phase factors e −iAn(Q) and e iAn(Q) into the nuclear wavefunctions ψ n (Q, t) and adiabatic states |n(Q) .
The time-dependent Schrödinger equation in the adiabatic representation, is obtained by substituting ansatz (2) into Eq.( 3) and projecting onto electronic states m(Q)| for m ∈ {1, . . ., S}.Note that we have introduced the representation-independent matrix notation: bold font indicates either an S × S matrix, i.e., an electronic operator, or S-dimensional vector, and the hat ˆindicates a nuclear operator.In particular, ( Ĥad ) mn = m|H|n is the adiabatic Hamiltonian, and ψ(t) denotes the molecular wavepacket in the adiabatic representation with components ψ n (t); henceforth, m, n ∈ {1, . . ., S} unless otherwise stated.The formal solution of Eq. ( 4) for a given initial condition ψ(0) is where Ûad (t) denotes the exact evolution operator Ûi (t) := exp(−i Ĥi t/ ) with i = ad.The adiabatic Hamiltonian Ĥad is often expressed as where we have used the matrix notation for the diagonal adiabatic potential energy matrix . The D-dimensional vector P is the nuclear momentum conjugate to Q and the dot • denotes a dot product in the D-dimensional nuclear vector space; we use the mass-scaled coordinates for simplicity.
Expressing the nonadiabatic vector couplings as shows that they are singular at a conical intersection, 65 which is a nuclear geometry Q 0 where 11,20,66 This singularity causes problems for the nonadiabatic dynamics simulations, and especially for the grid-based methods because an infinitely dense grid would be required to describe the singularity.The singularity, however, can be removed by a coordinate-dependent unitary transformation of Hamiltonian (6) into its quasidiabatic representation, where F qd (Q) and G qd (Q) denote the residual nonadiabatic vector and scalar couplings, respectively.Different transformation matrices S(Q) are obtained by different quasidiabatization schemes, 65 which include the block-diagonalization of the reference Hamiltonian matrix, 44,[67][68][69] integration of the nonadiabatic couplings, [70][71][72][73][74] use of the molecular properties, [75][76][77][78][79][80][81] and construction of regularized quasidiabatic states. 43,46,82We have chosen the regularized diabatization scheme because it is simple to implement and because it removes the conical intersection singularity reliably and efficiently. 83The choice of quasidiabatization affects the magnitude of the residual couplings and, therefore, their importance for the accuracy of nonadiabatic simulations. 40e initial state can be propagated with Ĥqd-exact instead of Ĥad to obtain the solution which is equivalent to ψ ad (t).However, it is much more common to use the simpler approximate quasidiabatic Hamiltonian This approximation is typically justified only heuristically by referring to the "small" magnitude of the residual nonadiabatic couplings.Nevertheless, the solution obtained using Ĥqd-approx is not exact and would still be only approximate even if evaluated numerically exactly.
We have numerically compared the three different solutions, ψ ad (t), ψ qd-exact (t), and ψ qd-approx (t) using the archetypal quadratic E ⊗ e Jahn-Teller model. 46,54In the model, two electronic states labeled by n = 1 and n = 2 are coupled by doubly degenerate normal modes Q 1 and Q 2 .We express the potential energy surface in polar coordinates-the radius in the space of the degenerate normal modes. 46,54Also, we work in natural units (n.u.) by setting k = M = = 1 n.u., where M is the mass associated with the degenerate normal modes, and ω = k/M = 1 n.u. is a quantum of the vibrational energy of these modes. 54,55e adiabatic potential energy surfaces are and Jahn-Teller coupling 46,54 The nonadiabatic vector coupling is 46,54 in our study, the coupling coefficients were c 1 = 1 n.u. and c 2 = 0.25 n.u.Unlike the potential energy V ad (Q), the nonadiabatic couplings F ad (Q) and G ad (Q) are affected by overall phases of the adiabatic states.One of the most standard choices holds exceptionally in the Jahn-Teller model and other systems in which a finite number of states represents the system exactly in both the adiabatic and diabatic representations.
Relationship (15) allows us to re-express Hamiltonian (6) in a simpler form which would generally only hold for S → ∞.
The exact quasidiabatic Hamiltonian is obtained from Hamiltonian (16) using the adiabatic to quasidiabatic transformation matrix 43,46,82 In Eq. ( 17), the (no longer diagonal) quasidiabatic potential energy matrix is and the residual nonadiabatic coupling is 42 where In realistic systems, Eqs. ( 15) and ( 16) would only hold if an infinite number of electronic states were considered.Likewise, the strictly diabatic Hamiltonian does not, in general, exist unless S → ∞. 41,42 Yet, due to its exceptional form, the Jahn-Teller Hamiltonian can be strictly diabatized 46 into a separable Hamiltonian obtained by replacing the transformation matrix S(Q) in Eq. ( 18) with the diabatic potential energy matrix is Potential energy surfaces V ad (Q), V qd (Q), and V diab (Q) in the vicinity of the conical intersection Q = 0 are visualized in Fig. 1.The strictly diabatic states are only coupled by the off-diagonal elements in V diab (Q) because the residual nonadiabatic couplings vanish: Like F ad (Q) and G ad (Q), the transformation matrices S(Q) and T(Q) in Eqs. ( 18) and ( 22) change according to overall phases of the adiabatic states (see Sec. S1 of the supplementary material).
We simulated the quantum dynamics following a transition from the ground vibrational eigenstate of the electronic state of A symmetry, FIG. 1. Potential energy surfaces in the E ⊗ e Jahn-Teller model around the conical intersection the two elements of the diagonal matrix V ad (Q) touch each other at the conical intersection.
where ψ(t) = T( Q)ψ(t) denotes the wavepacket in the strictly diabatic representation.To obtain the initial wavepacket, we assumed an impulsive excitation, i.e., the validity of the time-dependent perturbation theory and Condon approximation during the excitation.
Among various time propagation schemes, [33][34][35][36] we chose the geometric integrators 36,38,56 because they preserve exactly geometric properties of the exact evolution.Second-order split-operator algorithms, [37][38][39] including the TVT algorithm, preserve the linearity, norm, inner-product, symplecticity, stability, symmetry, and time reversibility of the exact solution. 36,38,56e implicit midpoint method, like the closely related trapezoidal rule (or Crank-Nicolson) method, 84,85 conserves, in addition, the energy.Both the split-operator and implicit midpoint methods can be symmetrically composed using various recursive or direct schemes 36,[57][58][59] to obtain integrators of arbitrary even orders of accuracy; 60,61 these compositions conserve all the geometric properties conserved by the elementary methods (see Refs. 60 and 61, and Sec.S2 of the supplementary material).The composed 36,57,58 TVT split-operator algorithm was used to propagate the wavepacket with the separable Hamiltonians ( Ĥdiab and Ĥqd-approx ), whereas the composed implicit midpoint method 56,84,85 was employed for propagations with the nonseparable Hamiltonians ( Ĥad and Ĥqd-exact ).Both integrators were composed using the optimal 59 eighth-order scheme, which, when combined with a small time step ∆t = 1/(40ω) = 0.025 n.u., led to time discretization errors negligible to the errors due to the use of different forms of the Hamiltonian (see Sec. S3 of the supplementary material).
On a grid of infinite range and density, nonadiabatic quantum dynamics simulated using Ĥdiab , Ĥad , and Ĥqd-exact would be identical.Therefore, the comparison of these Hamiltonians is only meaningful for a specific finite grid; we used a uniform grid of 64 × 64 points defined between Q l = −10 n.u. and Q l = 10 n.u. for l ∈ {1, 2}.Also, the favorable form of the strictly diabatic Hamiltonian allowed us to obtain the exact reference solution ψref (t) := Ûdiab (t) ψ(0) that is fully converged in both space and time: to ensure the grid convergence of the reference wavepacket, we used a grid of 128 × 128 points defined between In Sec.S3 of the supplementary material, we show that both the spatial and time discretization errors of ψref (t) are negligible (< 10 −10 ).In contrast, even on an infinite grid, the wavepacket ψqd-approx (t) would still not be exact because the residual nonadiabatic couplings are ignored.Section S3 of the supplementary material shows that even on a grid of 64×64 points, the spatial discretization errors are only minor contributors to the total errors of ψqd-approx (t).

III. RESULTS AND DISCUSSION
We compared the nonadiabatic quantum dynamics simulated using the adiabatic ( Ĥad ), exact quasidiabatic ( Ĥqd-exact ), and approximate quasidiabatic ( Ĥqd-approx ) Hamiltonians (see Sec. S4 of the supplementary material for the results obtained in the adiabatic representation without including the geometric phase).The reference quantum dynamics simulated using the strictly diabatic Hamiltonian Ĥdiab was left out of the comparison and was only used as the benchmark because Ĥdiab only exists, in general, when S → ∞.Before comparing the wavepackets themselves, we first present a comparison of three computed observables: the power spectrum obtained by Fourier transforming the autocorrelation function ψ(0)| ψ(t) (Fig. 2), population P 1 (t) of the first (n = 1) adiabatic electronic state (Fig. 3), and position ρ(t) (Fig. 4).The validity of this comparison is justified because the time discretization  .For applications that do not require extremely precise peak positions and intensities, even the spectrum obtained using Ĥqd-approx would suffice.In contrast, in panels (c) of Figs. 3 and 4, we see that both the population and position obtained with the approximate quasidiabatic Hamiltonian become inaccurate already after t ≈ 50 n.u.For t > 50 n.u., the population and position are described accurately only in simulations using either the exact quasidiabatic or strictly diabatic Hamiltonian.Among these Hamiltonians, however, only the exact quasidiabatic Hamiltonian exists in general unless S → ∞. 41,42 Some observables may be accurate, even if they are computed from a poor wavepacket.
In contrast, an accurate wavepacket ensures the accuracy of every observable computed from it.For a more stringent comparison between the different Hamiltonians, in Fig. 5 we, therefore, display the wavepackets ψ(t f ) at the final time.Whereas ψqd-exact (t) resembles  the exact wavepacket ψref (t) closely and ψqd-approx (t) has a similar overall shape but differs in the nodal structure and other details, ψad (t) is completely different.

IV. CONCLUSION
We rigorously compared the suitability of different forms of the molecular Hamiltonian for simulating the nonadiabatic quantum dynamics in the vicinity of a conical intersection.
This comparison was possible by taking advantage of the high-order geometric integrators and exceptional existence of the strictly diabatic Hamiltonian for the E ⊗ e Jahn-Teller model.The errors due to using the different forms of the molecular Hamiltonian were measured by comparing the fully converged exact reference simulation with simulations performed on a slightly sparser grid using the adiabatic, exact quasidiabatic, or approximate quasidiabatic Hamiltonian.We found that the nonadiabatic quantum dynamics simulated using the exact quasidiabatic Hamiltonian is nearly identical to the reference simulation obtained using the strictly diabatic Hamiltonian.The regularized diabatization scheme 43,46,82 was used for its simplicity, but the exact quasidiabatic Hamiltonian obtained through other schemes should lead to very similar and accurate results, as long as the quasidiabatization removes the conical intersection singularity, because Hamiltonian ( 8) is exact regardless of the quasidiabatization scheme.In contrast, the accuracy of the simulation with the approximate quasidiabatic Hamiltonian depends on the size of the neglected residual nonadiabatic couplings, and therefore, on the quasidiabatization scheme. 40 return to the question posed in the title, the approximate quasidiabatic Hamiltonian is appropriate if a quick solution of only moderate accuracy is required because the simple separable form of this Hamiltonian allows using more efficient time-propagation algorithms.
The accuracy can be further improved by employing more sophisticated quasidiabatization schemes, which reduce the size of the neglected residual couplings.The adiabatic Hamiltonian, although exact, is not suitable for simulating quantum dynamics at a conical intersection because of the singularity of the nonadiabatic couplings there; large errors appeared since this singularity could not be described well on a finite grid.8][89][90] Yet, our results clearly show that the rarely used exact quasidi-of these electronic states must be accompanied by compensating sign changes of the associated nuclear wavefunctions ψ dv n (Q, t) in order that the total molecular wavepacket |Ψ(Q, t) be single-valued. 22However, numerically representing, e.g., on a grid, a wavefunction that undergoes a sign change around a conical intersection is difficult.Another way of incorporating the geometric phase effect, according to the approach by Mead and Truhlar, 22 is to multiply the double-valued adiabatic electronic states by complex phase factors e iAn (Q) to obtain the equivalent adiabatic states |n(Q) = e iAn(Q) |n dv (Q) that are single-valued.
Following previous work, 25,26 we use the phase where k n must be an odd integer to compensate for the sign change of |n dv (Q) around a conical intersection.
The multiplication of |n dv (Q) by the coordinate-dependent phase factors does not affect the diagonal adiabatic potential energy matrix because e iAn(Q) commutes with the electronic Hamiltonian H e (Q): where The nuclear wavefunctions transform as ψ n (Q, t) = e −iAn(Q) ψ dv n (Q, t) and the nonadiabatic couplings as and In the Jahn-Teller model, the adiabatic potential energy matrix is defined by diagonalizing the strictly diabatic potential energy matrix of the model. 54In Refs.46 and 54, V diab (Q) is diagonalized by the transformation matrix which is one out of a continuum of possible transformation matrices that diagonalize This standard choice of the transformation matrix T dv (Q) yields the adiabatic electronic states that are coupled through the nonadiabatic vector couplings with vanishing diagonal elements.However, the states |n dv (Q) are double-valued because upon encircling a conical intersection [i.e., as the polar angle φ(Q) changes from −π to π], θ(Q) changes from −π/2 to π/2, causing |n dv (Q) to change sign.We, therefore, multiply the double-valued states |n dv (Q) with phase factors e iAn(Q) , where , to obtain the single-valued adiabatic electronic states These states are single-valued because θ + (Q) changes from −π to π and θ − (Q) starts and ends at zero as φ(Q) goes from −π to π.Moreover, like T dv (Q), the matrix T(Q) [see Eq. ( 22) of the main text], which transforms the strictly diabatic states |ñ directly into the single-valued adiabatic states |n(Q) , also diagonalizes V diab (Q) [into the same diagonal matrix as does T dv (Q)] because overall phases of the eigenvectors can be freely chosen (and here this phase can be a function of Q).The two transformation matrices T(Q) and T dv (Q) are related by ) δ mn is the diagonal matrix that transforms the double-valued states into the single-valued states, i.e., , where the matrix transforms the double-valued adiabatic states into the quasidiabatic states.

S2. CONSERVATION OF GEOMETRIC PROPERTIES BY THE TIME PROPAGATION SCHEMES
In Fig. S1, we show two of the geometric properties preserved exactly by the compositions of the implicit midpoint method: the norm of the wavepacket and expectation value of energy.Among the two properties, only the norm is conserved exactly by the compositions of the split-operator algorithm.See Refs.60 and 61 for complete analytical and numerical demonstrations of the conservation of linearity, inner-product, symplecticity, stability, symmetry, and time reversibility by symmetric compositions of the split-operator or implicit midpoint method.(h) i ∈ {qd-exact, qd-approx}, and E(t We used the optimal eighth-order composition 59 of the split-operator algorithm 37 to propagate the wavepacket with the separable Hamiltonians (either the approximate quasidiabatic or strictly diabatic Hamiltonian) and the optimal eighth-order composition 59 of the implicit midpoint method 56,85 to propagate the wavepacket with the nonseparable Hamiltonians (the adiabatic or exact quasidiabatic Hamiltonian).Because exact norm and energy conservation are built into the compositions of the implicit midpoint method, these geometric properties are conserved exactly regardless of the size of the time step.Likewise, the norm is conserved exactly by the compositions of the split-operator algorithm.The apparent conservation of energy by the split-operator algorithm in Fig. S1 only results from using a very small time step.In fact, it was shown in Ref. 61 that the energy conservation by the compositions of the split-operator algorithms is not exact but follows the order of convergence of the integrator.

S3. TIME AND SPATIAL DISCRETIZATION ERRORS
To quantify the errors due to different forms of the molecular Hamiltonian, an exact reference solution is required.In principle, the quantum dynamics simulated using any of the three exact Hamiltonians-the strictly diabatic, exact quasidiabatic, or adiabatic Hamiltonian-can serve as a benchmark, but the wavepacket propagated with the strictly diabatic Hamiltonian is the easiest to converge in both time and space because of the separable form of this Hamiltonian.We, therefore, chose the wavepacket propagated with the strictly diabatic Hamiltonian as the benchmark.In Fig. S2, we show that both the time and spatial discretization errors of propagation with the reference Hamiltonian are negligible (< 10 −10 ).In particular, the numerical errors of the reference wavepacket are minuscule compared to the errors due to different forms of the Hamiltonian (see Fig. 6 of the main text).We used the distance (time) ∆t [ψ(t)] := ψ (∆t,N ) (t) − ψ (∆t/2,N ) (t) to estimate the time discretization error of ψ (∆t,N ) (t); here, ψ (∆t,N ) (t) denotes the molecular wavepacket propagated to time t with the time step of ∆t on a grid of N × N points.The spatial discretization errors of ψ (∆t,N ) (t) were estimated by the distance The grid of 2N ×2N points is defined to be a factor of √ 2 wider in each dimension compared to the grid of N × N points.Correspondingly, the grid of 2N × 2N points is also a factor of √ 2 denser in each dimension compared to the grid of N × N points.For an observable A, we measured the time discretization errors using represented using the double-valued adiabatic states |n dv (Q) .However, like the electronic states themselves, the initial state ψ dv (0) = R( Q)ψ(0) in this representation is also doublevalued.We examine the effect of neglecting the geometric phase by simulating the dynamics with one of the two branches of the double-valued state ψ dv (0) as the initial wavepacket.
the quasidiabatic and strictly diabatic potential energy matrices are shown together because they are all equal.
errors of the presented observables and the time and spatial discretization errors of the reference observables are negligible (see Sec. S3 of the supplementary material).Panels (a) of Figs.2-4show that none of these observables is obtained accurately with the adiabatic Hamiltonian Ĥad even if the geometric phase is included: The positions and intensities of the peaks in the power spectrum are inaccurate, while the population P 1,ad (t) and position ρ(t) ad deviate very rapidly from their benchmark values.In contrast, all three presented observables are computed extremely accurately if the wavepacket is propagated with the exact quasidiabatic Hamiltonian [see panels (b) of Figs.2-4].There is no visible difference between the observables obtained using Ĥqd-exact and the benchmark observables obtained using Ĥdiab .

Figure 2 (
Figure 2(c) shows that the spectrum obtained by Fourier transforming ψ(0)| ψqd-approx (t) is very similar to the benchmark spectrum; the differences are only clearly visible in the zoomed-in version [see inset of Fig. 2(c)].For applications that do not require extremely

FIG. 3 .
FIG. 2. Power spectrum computed by Fourier transforming the autocorrelation function ψ(0)| ψi (t) of the wavepacket propagated with the (a) adiabatic (i = ad), (b) exact quasidiabatic (i = qd-exact), or (c) approximate quasidiabatic (i = qd-approx) Hamiltonian is compared with the benchmark spectrum (i = ref).To emulate the broadening of the peaks, the autocorrelation function was multiplied by the damping function f (t) = exp[−(t/t damp ) 2 ] with t damp = 80 n.u.before the Fourier transformation.Zoomed-in versions of the spectra are presented in the insets.

4 Q 1 FIG. 5 .FIG. 6 .
FIG. 5. Accuracy of the nonadiabatic quantum dynamics simulation demonstrated by comparing the wavepackets ψi (t f ) at the final time t = t f propagated with the (a) adiabatic (i = ad), (b) exact quasidiabatic (i = qd-exact), (c) approximate quasidiabatic (i = qd-approx), or (d) strictly diabatic (i = ref) Hamiltonian.The reference wavepacket ψref (t) serves as the benchmark.We only show Re[ ψ2,i (t f )], i.e., the real part of the nuclear wavepacket in the second (n = 2) electronic state of the strictly diabatic representation.
and the last step of Eq. (S1) holds since [V dv ad (Q)] mn = m dv (Q)|H e (Q)|n dv (Q) is diagonal.In contrast, both the nuclear wavefunctions and nonadiabatic couplings are affected by this transformation of the adiabatic electronic states:

FIG. S1 .
FIG. S1.Exact conservation of the norm ψ i (t) of the wavepacket [panels (a)-(d)] and of the expectation value E(t) i of energy [panels (e)-(h)].The initial value of the wavepacket norm is ψ i (0) = 1, for i ∈ {ad, qd-exact, qd-approx, ref}, and the initial values of the energy are E(0) ad = 2 n.u. and E(0) i = 1 n.u., for i ∈ {qd-exact, qd-approx, ref}.To evaluate the expectation value of energy, the wavepackets were transformed to the same representation as the Hamiltonian: E(t) ad := ψ ad (t)| Ĥad |ψ ad (t) , E(t) i := ψ i (t)|S( Q) † Ĥi S( Q)|ψ i (t) for FIG. S2.Time discretization errors (time) ∆t[ ψref (t)] (dotted line) and spatial discretization errors |A (∆t,N ) − A (∆t/2,N ) | and spatial discretization errors using(grid) N (A) := |A (∆t,N ) − A (∆t,2N ) |,where A (∆t,N ) is the observable obtained from a simulation on a grid of N × N points using the time step ∆t.FigureS3shows that the time and spatial discretization errors of all reference observables are negligible (< 3 • 10 −10 ) except for the spatial discretization errors of the population of the first (n = 1) adiabatic electronic state, which are almost entirely due to the spatial discretization errors of the unitary transformation of the wavepacket from the strictly diabatic to the adiabatic representation.Although larger, the spatial discretization errors (grid) N [P 1,ref (t)] of the reference population are still much smaller than the differences shown in Fig. 3 of the main text; even the virtually invisible differences between P 1,ref (t) and P 1,qd-exact (t) [in Fig. 3(b) of the main text] are approximately an order of magnitude greater than the spatial discretization errors of the reference population.

Figure
FigureS4shows that the time discretization errors of the propagation with the adiabatic, exact quasidiabatic, or approximate quasidiabatic Hamiltonian are negligible (i.e., at least two orders of magnitude smaller) compared to the errors due to different forms of the Hamiltonian (cf.Fig.6of the main text).For the time step ∆t = 0.025 n.u.employed in the simulations, the composed split-operator algorithm, used for propagations with the separable approximate quasidiabatic Hamiltonian, leads to much smaller time discretization FIG. S3.Time discretization errors (top) and the spatial discretization errors (bottom) of the reference observables: (a)-(b) the power spectrum, (c)-(d) population of the first (n = 1) adiabatic electronic state, and (e)-(f) position.Line labels are the same as in Fig. S2.
FIG. S4.Time discretization errors of the wavepacket propagated with the (a) adiabatic (i = ad), (b) exact quasidiabatic (i = qd-exact), or (c) approximate quasidiabatic (i = qd-approx) Hamiltonian.The eighth-order composition of the implicit midpoint method was used to propagate the wavepacket with either the adiabatic or exact quasidiabatic Hamiltonian.The eighth-order composition of the TVT split-operator algorithm was used to propagate the wavepacket with the approximate quasidiabatic Hamiltonian.The time step ∆t = 0.025 n.u. was used in all simulations.

Figure
Figure S7 shows that the power spectrum I i (ω) [panel (a)], population P 1,i (ω) [panel (b)], and position ρ(t) i [panel (c)] obtained in the adiabatic representation without the geometric phase (i = ad-no-gp) are much less accurate than those obtained in the adiabatic representation in which the geometric phase is included [i = ad, Figs.2(a), 3(a), and 4(a) of the main text].