How important are the residual nonadiabatic couplings for an accurate simulation of nonadiabatic quantum dynamics in a quasidiabatic representation?

The celebrated Born–Oppenheimer approximation,¹ which treats the electronic and nuclear motions in molecules separately, is no longer valid for describing processes involving two or more strongly vibronically coupled electronic states. A common approach that goes beyond this approximation^2–9 consists in solving the time-dependent Schrödinger equation with a truncated molecular Hamiltonian that includes only a few most significantly coupled^10,11 Born–Oppenheimer electronic states.^12–14 The “adiabatic” states, obtained directly from the electronic structure calculations, are, however, not adequate for representing the molecular Hamiltonian in the region of strong nonadiabatic couplings; in particular, the couplings between the states diverge at conical intersections,^2,14–19 where potential energy surfaces of two or more adiabatic states intersect.

Quasidiabatization, i.e., a coordinate-dependent unitary transformation^20–22 of the molecular Hamiltonian that reduces the magnitude of the nonadiabatic vector couplings, rectifies this singularity. The transformation matrix can be obtained by various quasidiabatization schemes, of which a few representative examples include methods based on the integration of the nonadiabatic couplings^23–33 or on different molecular properties^34–40 and the block-diagonalization^41–44 or regularized diabatization^45–47 schemes.

In systems with more than one nuclear degree of freedom, the strict diabatization, which eliminates the nonadiabatic couplings completely, is only possible if infinitely many electronic states are considered.^21,48 The best one can do for a general subsystem with a finite number of electronic states is the above-mentioned quasidiabatization in which the unitary transformation reduces the magnitude of the couplings but does not remove them entirely. However, it is a common practice to neglect these nonvanishing “residual” couplings present in the exact quasidiabatic Hamiltonian and thus obtain an approximate quasidiabatic Hamiltonian, whose additional benefit is a simpler, separable form convenient for quantum simulations.

Here, we propose a general method that quantifies the importance of the residual couplings by comparing nonadiabatic simulations performed either with the exact quasidiabatic Hamiltonian—obtained through an exact unitary transformation of the adiabatic Hamiltonian—or with the approximate quasidiabatic Hamiltonian, which neglects the residual couplings. By definition and regardless of the magnitude of the residual couplings, the results obtained with the exact quasidiabatic Hamiltonian can serve as the exact benchmark, as long as the numerical errors are negligible.⁴⁹ Therefore, for a valid comparison, one needs a time propagation scheme that can treat even the nonseparable exact quasidiabatic Hamiltonian and that ensures that the numerical errors are negligible to the errors due to neglecting the residual couplings. Among various integrators^50–57 that satisfy both requirements, we chose the optimal eighth-order⁵⁸ composition^57,59–61 of the implicit midpoint method^56,57,62 because it also preserves exactly⁵⁴ various geometric properties of the exact solution.^56,57

After presenting the general method in Sec. II, in Sec. III, we provide realistic numerical examples in which we employ the method to quantify the importance of the residual couplings in nonadiabatic simulations of nitrogen trioxide (NO₃)^63–66 and hydrogen cyanide (HCN).^41,67–69 Whereas the NO₃ model was quasidiabatized with the regularized diabatization scheme,^45–47 the block-diagonalization scheme^41–44 was employed in the HCN model. To find out how the errors due to ignoring the residual couplings depend on the sophistication of the quasidiabatization and on the initial state, in Sec. III C, we compare the first- and second-order regularized diabatization schemes^45–47 on the model of a displaced excitation of NO₃.

We begin by introducing the standard molecular Hamiltonian $H=TN+Te+V$⁠, where $TN$ and $Te$ are the kinetic energy operators of the nuclei and electrons and $V$ is the molecular potential energy operator. One may express the molecular Hamiltonian equivalently as $H=TN+He$ by defining the electronic Hamiltonian $He≔Te+V$⁠, an operator acting on the electronic degrees of freedom and depending parametrically on the nuclear coordinates, described by a D-dimensional vector Q. For each fixed nuclear geometry, the time-independent Schrödinger equation,

for $He(Q)$ can be solved to obtain the nth adiabatic electronic state |n(Q)⟩ and potential energy surface V_n(Q) for $n∈N$⁠.

The adiabatic electronic eigenstates |n(Q)⟩, which depend on the nuclear coordinates Q, form a complete orthonormal set and can be employed to expand the exact solution of the time-dependent molecular Schrödinger equation,

with Hamiltonian $H$ as an infinite series

Note that Eqs. (2) and (3) combine the coordinate representation for the nuclei with the representation-independent Dirac notation for the electronic states; $ψnad(Q,t)$ is the time-dependent nuclear wavefunction (a wavepacket) on the nth adiabatic electronic surface. The Born–Huang expansion⁷⁰ of Eq. (3) is exact when an infinite number of electronic states are included, but in practice, |Ψ(Q,t)⟩_exact is approximated by truncating the sum in Eq. (3) and including only the most important S electronic states:^10,11,13

For brevity, we shall omit the subscript “trunc” in |Ψ(Q,t)⟩_trunc from now on.

Substituting ansatz (4) into the time-dependent Schrödinger equation (2) and projecting onto states ⟨m(Q)| for m ∈ {1, …, S} leads to the ordinary differential equation

expressed in a compact representation-independent matrix notation: bold font indicates either an S × S matrix (i.e., an electronic operator) or an S-dimensional vector and the hat (⁠$^$⁠) denotes a nuclear operator. In particular, $H^ad$ is the adiabatic Hamiltonian matrix with elements $(H^ad)mn=〈m|H|n〉$ and ψ^ad(t) is the molecular wavepacket in the adiabatic representation with components $ψnad(t)$⁠. Assuming the standard form $TN=P^2/2M$ of the nuclear kinetic energy operator, the adiabatic Hamiltonian matrix is given by the formula^2,4,14,21,71

which depends on the diagonal adiabatic potential energy matrix $[Vad(Q)]mn≔Vn(Q)δmn$⁠, the nonadiabatic vector couplings $[Fad(Q)]mn≔〈m(Q)|∇n(Q)〉$⁠, and the nonadiabatic scalar couplings $[Gad(Q)]mn≔〈m(Q)|∇2n(Q)〉$⁠. The dot (·) denotes the dot product in the D-dimensional nuclear vector space, and P is the canonical momentum conjugate to Q. Note that, for simplicity, the nuclear coordinates have been scaled so that each nuclear degree of freedom has the same mass M and, therefore, M is a scalar.

In practice, the nonadiabatic scalar couplings G_ad(Q) in Eq. (6) are often neglected, but this approximation can cause significant errors;^72,73 for the adiabatic Hamiltonian to be exact, it must include both F_ad(Q) and G_ad(Q). In Eqs. (3)–(6), one can freely choose overall phases of the adiabatic electronic states because both |n(Q)⟩ and $eiAn(Q)|n(Q)〉$ [where A_n(Q) are coordinate-dependent phases] are orthonormalized solutions of Eq. (1). In Ref. 49, we show how the choice of A_n(Q) affects the nonadiabatic couplings F_ad(Q) and G_ad(Q); in contrast, V_ad(Q) remains unaffected.

The nonadiabatic vector couplings can be re-expressed using the Hellmann–Feynman theorem as

accentuating the singularity of these couplings at a conical intersection^3,20—a nuclear geometry Q₀ where V_m(Q₀) = V_n(Q₀) for m ≠ n.^2,14,19 Moreover, Meek and Levine⁷⁴ pointed out that, unlike the singularity of $[Fad(Q)]mn$⁠, the singularity in the diagonal elements $[Gad(Q)]nn$ of the nonadiabatic scalar couplings is not even integrable over domains containing a conical intersection. Another complication associated with conical intersections is the geometric phase effect: the sign change of the real-valued adiabatic electronic state |n(Q)⟩ when transported along a loop containing a conical intersection.^75–87 Although the geometric phase effect can be effectively incorporated into nonadiabatic simulations in the adiabatic basis by appropriately choosing the above-mentioned phases A_n(Q) so that the states are single-valued,^75–87 the numerically problematic singularity remains.⁴⁹ Yet, both complications, namely, the singularity of the nonadiabatic couplings and the geometric phase effect, can be avoided simultaneously by transforming the adiabatic Hamiltonian to the quasidiabatic basis,

and thus obtaining the exact quasidiabatic Hamiltonian,

where $[Vqd(Q)]mn≔〈m′(Q)|He(Q)|n′(Q)〉$ is the nondiagonal quasidiabatic potential energy matrix, while $[Fqd(Q)]mn≔〈m′(Q)|∇n′(Q)〉$ and $[Gqd(Q)]mn≔〈m′(Q)|∇2n′(Q)〉$ are the residual vector and scalar couplings, respectively. Note that the molecular state |Ψ(Q, t)⟩ from Eq. (4) is independent of the choice of basis because the transformation (8) from the adiabatic to quasidiabatic electronic basis is accompanied by a simultaneous transformation,

of nuclear wavefunctions. The transformation matrix S(Q) is obtained by any of the many quasidiabatization schemes,^20–47 but the magnitude of the residual nonadiabatic couplings depends on the scheme. Following Ref. 21, we measure this magnitude with the quantity

where

is the square of the Frobenius norm of F_qd(Q) [note that the evaluation of Eq. (12) involves both the matrix product of S × S matrices and the scalar product of D-vectors]. Section S1 of the supplementary material describes the numerical evaluation of $R[Fqd(Q)]$ in further detail.

It is well-known that, unless S is infinite or D = 1, in a general system, no diabatization scheme yields the strictly diabatic states [i.e., states in which the exact Hamiltonian (9) has zero residual nonadiabatic couplings].^21,48 The transformation by a finite S × S matrix S(Q) can only lead to quasidiabatic states, which are coupled both by the off-diagonal (m ≠ n) elements $[Vqd(Q)]mn$ of the quasidiabatic potential energy matrix and by the—perhaps small but nonvanishing—residual nonadiabatic couplings.⁸⁸ In practice, however, these residual couplings are often ignored in Eq. (9) in order to obtain the approximate quasidiabatic Hamiltonian,

Although the magnitude $R[Fqd(Q)]$ itself may indicate whether it is admissible to neglect the residual couplings, a much more rigorous way to quantify the impact of this approximation on a particular nonadiabatic simulation consists in evaluating the quantum fidelity,⁸⁹ 

and distance,

between the states ψ_qd-approx(t) and ψ_qd-exact(t), evolved with the approximate and exact quasidiabatic Hamiltonians, respectively [i.e., $ψi(t)=exp(−iH^it/ℏ)ψ(0)$ for i ∈ {qd-approx, qd-exact}]. The more important the residual couplings, the smaller the quantum fidelity and the larger the distance.

By both propagating and comparing the wavepackets ψ_qd-approx(t) and ψ_qd-exact(t) in the same quasidiabatic representation, one avoids contaminating the errors due to the neglect of the residual couplings with the numerical errors due to the transformation between representations. In fact, as long as it is numerically converged, ψ_qd-exact(t) serves as the exact benchmark regardless of the size of the residual couplings because the exact quasidiabatic and adiabatic Hamiltonians are exact unitary transformations of each other.⁴⁹ 

The exact quasidiabatic Hamiltonian from Eq. (9) cannot be expressed as a sum of terms depending purely on either the position or the momentum operator. Due to this nonseparable nature of the Hamiltonian, we require an integrator that is applicable to any form of the Hamiltonian. For example, the popular split-operator algorithm^61,90–92 cannot be employed. The wavepackets are, therefore, propagated with the composition^57,59–61 of the implicit midpoint method,^56,57,62 which, like the closely related trapezoidal rule (or Crank–Nicolson method),^62,93 works for both separable and nonseparable Hamiltonians, as long as the action of the Hamiltonian on the wavepacket can be evaluated. Moreover, in contrast to some other methods applicable to nonseparable Hamiltonians, the chosen methods preserve exactly most geometric properties of the exact solution: norm, energy, and inner-product conservation, linearity, symplecticity, stability, symmetry, and time reversibility.⁵⁴ 

For a valid comparison of the two wavepackets propagated with either the exact or approximate quasidiabatic Hamiltonian, the numerical errors must be much smaller than the errors due to omitting the residual couplings. Owing to its exact symmetry, the implicit midpoint method can be composed using various schemes^57–60 to obtain integrators of arbitrary even orders of accuracy in the time step;^54,92 we compose the implicit midpoint method according to the optimal scheme⁵⁸ to obtain an eighth-order integrator. By using this high-order integrator with a small time step, the time discretization errors are kept negligible (see Sec. S2 of the supplementary material).

We now apply the method proposed in Sec. II to nonadiabatic quantum simulations in the cubic E ⊗ e Jahn–Teller model of NO₃^63–66 and in the induced Renner–Teller model of HCN.^41,67–69 Despite their reduced dimensionality, these two-dimensional (D = 2) two-state (S = 2) models exhibit interesting dynamics^64,69,94 due to the presence of strong nonadiabatic couplings; in particular, F_ad(Q) diverges at Q = 0, the point of intersection between the two adiabatic potential energy surfaces.

In both models, doubly degenerate electronic states labeled n = 1 and n = 2 are coupled by doubly degenerate normal modes Q₁ and Q₂. We use “natural” units (n.u.) throughout by setting k = M = ℏ = 1 n.u., where M is the mass associated with the degenerate normal modes (which differs from the electron mass used in atomic units) and $ℏω=ℏk/M$ is the quantum of the vibrational energy of these modes. Whenever convenient, we express the potential energy surface in polar coordinates: the radius $ρ(Q)≔Q12+Q22$ and polar angle ϕ(Q) ≔ arctan(Q₂/Q₁).

All numerical wavepacket propagations were performed with a small time step of Δt = 1/(40ω) = 0.025 n.u. on a uniform grid of N × N points defined between Q_l = −Q_lim and Q_l = Q_lim in both nuclear dimensions: N = 64 and Q_lim = 10 n.u. in the NO₃ model, while N = 32 and Q_lim = 7 n.u. in the HCN model.

Although the strictly diabatic Hamiltonian,

does not exist in general, it may exist exceptionally and, in fact, is used to define the Jahn–Teller model.^45,63–66 In Eq. (16), the diabatic potential energy matrix,

depends on the cubic potential energy E₀(ρ, ϕ) ≔ kρ²/2 + 2αρ³ cos 3ϕ and Jahn–Teller coupling,⁶⁴ 

where f(ρ) ≔ c₁ρ + c₃ρ³. In our nonadiabatic simulations of nitrogen trioxide, we used the Jahn–Teller model of NO₃ from Ref. 64 with parameters α = −0.0125 n.u., c₁ = 0.375 n.u., c₂ = −0.0668 n.u., and c₃ = −0.0119 n.u. To simplify the following presentation, we rewrite E_cpl(Q) as E_cpl(Q) = |E_cpl(Q)|e^−2iθ(Q) using the mixing angle

Our previous study⁴⁹ on a similar system showed that the exact quasidiabatic and strictly diabatic Hamiltonians yield nearly identical results. Here, however, we intentionally avoid using the strictly diabatic Hamiltonian as a benchmark and use it only to define the model, in order that the approach and conclusions of this study are applicable also to systems where the strictly diabatic Hamiltonian does not exist^21,48 (see Sec. III B for an explicit example of such a system).

The adiabatic states in the Jahn–Teller model are obtained by a process inverse to diabatization, i.e., by a unitary transformation of the strictly diabatic states using any matrix that diagonalizes V_diab(Q). Following Refs. 45 and 63, we employed the transformation matrix

In the resulting adiabatic representation, the diagonal potential energy matrix has elements V₁(Q) = V₊(Q) and V₂(Q) = V₋(Q), where V_±(Q) ≔ E₀(Q) ± |E_cpl(Q)|. [Matrix elements of V_ad(Q) and V_diab(Q) are plotted in Fig. 1]. Transformation (20) also yields analytical expressions for the nonadiabatic vector couplings, ^63,64

and for the nonadiabatic scalar couplings,

As expected, the nonadiabatic couplings diverge at the conical intersection at ρ = 0 since the azimuthal component of F_ad(Q) is proportional to

In the cubic Jahn–Teller model, the regularized diabatization scheme^45–47 can be implemented analytically. The jth order adiabatic to quasidiabatic transformation matrix,

is obtained simply by replacing θ(Q) with θ^(j)(Q) in Eq. (20) for T(Q): θ⁽¹⁾(ρ, ϕ) ≔ ϕ/2 in the first-order scheme and

in the second-order scheme, while—in the cubic Jahn–Teller model—the third-order quasidiabatization is already identical to the strict diabatization, i.e., θ⁽³⁾(ρ, ϕ) = θ(ρ, ϕ). The quasidiabatization yields the potential energy matrix,

residual vector couplings,

and residual scalar couplings,

where $θ−(j)(Q)≔θ(Q)−θ(j)(Q)$⁠. The resulting magnitude of the residual couplings in the first-order (j = 1) scheme is $R[Fqd(Q)]=3.8$ n.u.

The hermiticity of Hamiltonian (9) is broken on a finite grid because the commutator relation $[P^,Fqd(Q^)]=−iℏ∇⋅Fqd(Q^)$ holds only approximately unless the grid is infinitely dense. Yet, the hermiticity of the Hamiltonian is essential for the norm conservation (see Fig. S5 in Sec. S3 of the supplementary material), which, in turn, is required for quantum fidelity $F(t)$ and distance $D(t)$ to be valid measures of the importance of the residual nonadiabatic couplings. To make the exact quasidiabatic Hamiltonian exactly Hermitian, we re-express it as

using the relationship

which holds—exceptionally—for systems, such as the Jahn–Teller model, that can be represented exactly by a finite number of states; in general, Eq. (30) only holds when S → ∞.

Another benefit of Hamiltonian (29) is the absence of G_qd(Q), the evaluation of which represents the computational bottleneck in realistic systems. Likewise, the equations of motion in the widely employed Meyer–Miller approach^95–98 can be simplified greatly⁷² by starting from Hamiltonian (29) instead of Hamiltonian (9). These two forms of the molecular Hamiltonian, however, are strictly equivalent only if the electronic basis is complete.^21,71 In generic systems, in which the relation (30) does not hold and one is obliged to use the original Hamiltonian (9) and evaluate G_qd(Q), the computationally expensive evaluation of the second derivatives of electronic wavefunctions with respect to nuclear coordinates can still be avoided by using the relation

where $[Kqd(Q)]mn≔〈∇m′(Q)|∇n′(Q)〉$ requires only the first derivatives of the quasidiabatic electronic states |n′(Q)⟩ introduced in Eq. (8). In contrast to Eq. (30), relation (31) holds in arbitrary systems and for finite S.

To analyze the importance of residual couplings in NO₃, we simulated, with either the exact or approximate quasidiabatic Hamiltonian, the quantum dynamics following an electronic transition from the ground vibrational eigenstate of the ground electronic state V_g(Q) = −E_gap + kρ(Q)²/2 with E_gap = 11 n.u. (1 n.u. of energy here corresponds to 0.2 eV ≈ 0.007 a.u.). Invoking the time-dependent perturbation theory and Condon approximation, we considered the initial state in the quasidiabatic representation to be

where we omitted (also will omit) the superscript “qd” on the wavepacket for brevity.

Figure 2 shows that, in the nonadiabatic dynamics following the vertical excitation of NO₃, neglecting the residual couplings does not significantly affect the wavepacket [compare panels (a) and (b)], power spectrum I(ω) [panel (c)], or population $P1ad(t)$ [panel (d)]. Even the fidelity $F(t)$ [panel (e)] between the wavepackets propagated either with or without the residual couplings remains close to the maximal value of 1 until the final time t_f. Section S4 of the supplementary material further supports this conclusion by displaying the time dependence of position ⟨ρ⟩(t), potential energy ⟨V_qd⟩(t), and distance $D(t)$⁠. In contrast, as we will see in Secs. III B and III C, the residual couplings are much more significant in the HCN model and in the displaced excitation of NO₃.

In Sec. S5 of the supplementary material, we also analyze the importance of the residual couplings for different Jahn–Teller coupling coefficients and different initial populations.

The model of the induced Renner–Teller effect^41,67–69 is more realistic than the Jahn–Teller model from Sec. III A: In particular, the strictly diabatic Hamiltonian (16) cannot be defined and relationship (30) does not hold. Nevertheless, similar to the Jahn–Teller model, the nonadiabatic couplings between the adiabatic states are singular at Q = 0.⁴¹ Since Eq. (30) does not hold in the induced Renner–Teller model, the exactly Hermitian Hamiltonian (29) cannot be used instead of Hamiltonian (9). Yet, even with Hamiltonian (9), the norm is sufficiently converged in the grid density for the quantum fidelity $F(t)$ and distance $D(t)$ to be valid (see Fig. S6 of the supplementary material).

We follow Ref. 41, where the induced Renner–Teller model is quasidiabatized with the block-diagonalization scheme, which minimizes the residual couplings locally (around Q = 0 in this model).²¹ The resulting quasidiabatic potential energy matrix,

with V_±(ρ) ≔ [V₁(ρ) ± V₂(ρ)]/2 depends on the adiabatic potential energy surfaces V₁(ρ) = Δ + E_h(ρ) − $w$(ρ) and V₂(ρ) = E_h(ρ), where E_h(ρ) ≔ kρ²/2 and $w(ρ)≔(Δ2+2λ2ρ2)1/2$⁠. Analytical expressions for the nonadiabatic couplings, F_ad(Q) and G_ad(Q), and adiabatic to quasidiabatic transformation matrix S(Q) can be found in Ref. 41. In our nonadiabatic simulations of hydrogen cyanide, we used the induced Renner–Teller model of HCN from Refs. 41 and 67 with parameters Δ = 1.11 n.u. and λ = 1 n.u. The residual vector and scalar couplings⁴¹ are

and

respectively, where we have defined the D-dimensional (here, D = 2) vector $v$_F(Q) with components $v$_F(Q)₁ = −f_F(ρ(Q))Q₂ and $v$_F(Q)₂ = f_F(ρ(Q))Q₁ and functions f_F(ρ) ≔ [1 − $w$₊(ρ)]/ρ² and

with $w±(ρ)≔[1±Δ/w(ρ)]/2$⁠. The magnitude of the residual couplings (34) is $R[Fqd(Q)]=0.37$ n.u.; the adiabatic potential energy surfaces and functions f_F and f_G are plotted in Fig. 3.

Similar to Sec. III A, we simulate the dynamics following an electronic transition from the ground vibrational eigenstate of V_g(ρ, ϕ) = −E_gap + E_h(ρ) with E_gap = 153 n.u. (1 n.u. of energy here corresponds to 0.09 eV ≈ 0.003 a.u.). Unlike their analogs in Sec. III A, however, the two wavepackets, propagated with either the exact or approximate quasidiabatic Hamiltonian, differ significantly [compare panels (a) and (b) of Fig. 4]. Ignoring the residual couplings also leads to large errors in the power spectrum I(ω) [panel (c)], population $P1ad(t)$ [panel (d)], and fast decay of quantum fidelity $F(t)$ [panel (e)]. In particular, the population obtained with the approximate quasidiabatic Hamiltonian cannot be trusted because, e.g., at t = 169 n.u., the error $ϵres-cpl[P1ad(t)]≔|P1,qd-approxad(t)−P1,qd-exactad(t)|$ due to the neglect of the residual couplings is of the same order as the range $RP1ad≔P1,maxad−P1,minad$ of the population in the whole simulation interval: $ϵres-cpl[P1ad(t)]/RP1ad=0.7$⁠. Note also that neither $P1ad(t)$ nor $F(t)$ is affected by the overall phases of the two wavepackets, although a linearly growing overall phase difference appears to be the main contribution to the change of the spectrum, which is mostly shifted [see panel (c) of Fig. 4].

Although the excited states of NO₃ from Sec. III A are bright states, the coupled states modeled by the cubic E ⊗ e Jahn–Teller Hamiltonian can sometimes be dark. A wavepacket might reach such dark states at a nuclear geometry that is not the ground state equilibrium (e.g., via an intersection with a bright state). Motivated by this observation from Ref. 64, we consider the initial state (32) displaced in both Q₁ and Q₂ by −0.8 n.u. Keeping all other parameters fixed as in Sec. III A will allow us to analyze how the importance of the residual couplings depends on the initial state and on the quasidiabatization method used.

Figure 5 shows that, in contrast to the vertical excitation from Sec. III A (analyzed in Fig. 2), ignoring the residual couplings obtained by applying the first-order regularized scheme to the displaced excitation of NO₃ affects both the wavepacket [compare panels (a) and (b)] and observables significantly. Neither the spectrum I(ω) [panel (c)] nor population $P1ad(t)$ [panel (d)] is obtained accurately with the approximate quasidiabatic Hamiltonian. For example, at t = 25 n.u., the error of the population is almost half of the range of the population in the whole simulation interval: $ϵres-cpl[P1ad(t)]/RP1ad=0.4$⁠. The quantum fidelity [panel (e)] decreases rapidly to $F(tf)≈0.3$ at t_f = 50 n.u.

The residual couplings, however, can be made less important by an improved quasidiabatization. One can reduce the magnitude of the residual couplings from $R[Fqd(1)(Q)]=3.8$ n.u. to $R[Fqd(2)(Q)]=0.5$ n.u. by employing the more sophisticated second-order regularized diabatization scheme^45–47 obtained by inserting θ⁽²⁾(Q) from Eq. (25) into Eqs. (24) and (26)–(28). When this second-order scheme is used, the errors of the wavepacket ψ(t), spectrum I(ω), and population $P1ad(t)$ due to the neglect of the residual couplings all remain small (see Fig. 6); in particular, quantum fidelity $F(t)$ remains above 0.95 for all times until the final time t_f = 50 n.u. [see panel (e)]. {Note that the exact benchmark wavepackets [in panels (a) of Figs. 5 and 6] propagated in the two different quasidiabatic representations are slightly different not only because they are displayed in different representations but also because the initial states are different—they have the same analytical form but in two different quasidiabatic representations.}

We have shown that the common practice of neglecting the residual nonadiabatic couplings between quasidiabatic states can substantially lower the accuracy of nonadiabatic simulations and that the decrease in accuracy depends on the system, initial state, and employed quasidiabatization scheme. One can, therefore, answer the question posed in the title only after a careful analysis. In Sec. III, we have provided several examples where the approximate quasidiabatic Hamiltonian gives wrong results. Because it is potentially dangerous to employ an approximation without evaluating its impact, we have proposed a method to rigorously quantify the errors caused by ignoring the residual couplings.

When the residual couplings are significant and cannot be neglected, we suggest performing nonadiabatic simulations with the rarely used exact quasidiabatic Hamiltonian (9), which not only is analytically equivalent to the adiabatic Hamiltonian (6) but also yields numerically accurate results regardless of the magnitude of the residual couplings (as shown in Sec. S2 of the supplementary material and in Ref. 49). Although the general applicability of the exact quasidiabatic Hamiltonian depends on the availability of residual nonadiabatic couplings, these can be evaluated by employing recently developed schemes^99–103 even in rather complicated multi-state systems involving multiple conical intersections (including those between three electronic states^104–108). In complex systems where all practical quasidiabatization schemes lead to significant residual couplings, propagating the wavepacket with the exact quasidiabatic Hamiltonian would be particularly beneficial. Although the nonseparable form of this Hamiltonian complicates the time propagation, there exist efficient geometric integrators, such as the high-order compositions of the implicit midpoint method used here, which are applicable even to such Hamiltonians.

Last but not least, an accurate propagation of the wavepacket with the exact quasidiabatic Hamiltonian would be extremely useful for establishing highly accurate benchmarks in unfamiliar systems, where the impact of the residual nonadiabatic couplings on the quantum dynamics simulations is not yet known.

See the supplementary material for the details of the numerical evaluation of the magnitude of the residual couplings (Sec. S1), demonstration of the negligibility of spatial and time discretization errors (Sec. S2), conservation of geometric properties by the implicit midpoint method (Sec. S3), time dependence of position, potential energy, and distance (Sec. S4), and importance of the residual couplings for different Jahn–Teller coupling coefficients and different initial populations (Sec. S5).

The authors acknowledge the financial support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant Agreement No. 683069-MOLEQULE) and thank Tomislav Begušić and Nikolay Golubev for useful discussions.

The data that support the findings of this study are available within the article and its supplementary material.