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

Diabatization of the molecular Hamiltonian is a standard approach to remove the singularities of nonadiabatic couplings at conical intersections of adiabatic potential energy surfaces. In general, it is impossible to eliminate the nonadiabatic couplings entirely—the resulting “quasidiabatic” states are still coupled by smaller but nonvanishing residual nonadiabatic couplings, which are typically neglected. Here, we propose a general method for assessing the validity of this potentially drastic approximation by comparing quantum dynamics simulated either with or without the residual couplings. To make the numerical errors negligible to the errors due to neglecting the residual couplings, we use the highly accurate and general eighth-order composition of the implicit midpoint method. The usefulness of the proposed method is demonstrated on nonadiabatic simulations in the cubic Jahn–Teller model of nitrogen trioxide and in the induced Renner–Teller model of hydrogen cyanide. We find that, depending on the system, initial state, and employed quasidiabatization scheme, neglecting the residual couplings can result in wrong dynamics. In contrast, simulations with the exact quasidiabatic Hamiltonian, which contains the residual couplings, always yield accurate results. © 2021 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/). https://doi.org/10.1063/5.0046067., s


I. INTRODUCTION
The celebrated Born-Oppenheimer approximation, 1 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][3][4][5][6][7][8][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][13][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][15][16][17][18][19] where potential energy surfaces of two or more adiabatic states intersect.
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 Hamiltonianobtained 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. 49 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][51][52][53][54][55][56][57] that satisfy both requirements, we chose the optimal eighth-order 58 composition 57,59-61 of the implicit midpoint method 56,57,62 because it also preserves exactly 54 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 3 ) 63-66 and hydrogen cyanide (HCN). 41,[67][68][69] Whereas the NO 3 model was quasidiabatized with the regularized diabatization scheme, 45-47 the blockdiagonalization scheme [41][42][43][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 3 .

II. THEORY
We begin by introducing the standard molecular Hamiltonian H = T N + Te + V, where T N 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 = T N + 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 Vn(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; ψ ad n (Q, t) is the time-dependent nuclear wavefunction (a wavepacket) on the nth adiabatic electronic surface. The Born-Huang expansion 70 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,Ĥ ad is the adiabatic Hamiltonian matrix with elements (Ĥ ad )mn = ⟨m|H|n⟩ and ψ ad (t) is the molecular wavepacket in the adiabatic representation with components ψ ad n (t). Assuming the standard form T N =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 [V ad (Q)]mn ∶= Vn(Q)δmn, the nonadiabatic vector couplings [F ad (Q)]mn ∶= ⟨m(Q)|∇n(Q)⟩, and the nonadiabatic scalar couplings [G ad (Q)]mn ∶= ⟨m(Q)|∇ 2 n(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. 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 0 where Vm(Q 0 ) = Vn(Q 0 ) for m ≠ n. 2,14,19 Moreover, Meek and Levine 74 pointed out that, unlike the singularity of [F ad (Q)]mn, the singularity in the diagonal elements [G ad (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][76][77][78][79][80][81][82][83][84][85][86][87] Although the geometric phase effect can be effectively incorporated into nonadiabatic simulations in the adiabatic basis by appropriately choosing the above-mentioned phases An(Q) so that the states are single-valued, 75-87 the numerically problematic singularity remains. 49 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, ⟩ 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 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[F qd (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 [V qd (Q)]mn of the quasidiabatic potential energy matrix and by the-perhaps small but nonvanishing-residual nonadiabatic couplings. 88 In practice, however, these residual couplings are often ignored in Eq. (9) in order to obtain the approximate quasidiabatic Hamiltonian, Although the magnitude R[F qd (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, 89 and distance, between the states ψ qd-approx (t) and ψ qd-exact (t), evolved with the approximate and exact quasidiabatic Hamiltonians, respectively 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. 49 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. 54 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 The Journal of Chemical Physics ARTICLE scitation.org/journal/jcp 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 58 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).

III. NUMERICAL EXAMPLES
We now apply the method proposed in Sec. II to nonadiabatic quantum simulations in the cubic E ⊗ e Jahn-Teller model of NO 3 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 1 and Q 2 . We use "natural" units (n.u.) throughout by setting k = M = ̵ h = 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 ̵ hω = ̵ h √ 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) ∶= √ Q 2 1 + Q 2 2 and polar angle ϕ(Q) ∶= arctan(Q 2 /Q 1 ). 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 3 model, while N = 32 and Q lim = 7 n.u. in the HCN model.

A. Jahn-Teller effect in nitrogen trioxide
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][64][65][66] In Eq. (16), the diabatic potential energy matrix, depends on the cubic potential energy E 0 (ρ, ϕ) ∶= kρ 2 /2 + 2αρ 3 cos 3ϕ and Jahn-Teller coupling, 64 where f (ρ) ∶= c 1 ρ + c 3 ρ 3 . In our nonadiabatic simulations of nitrogen trioxide, we used the Jahn-Teller model of NO 3 Our previous study 49 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 hermiticity of Hamiltonian (9) is broken on a finite grid because the commutator relation [P, F qd (Q)] = −i ̵ h∇ ⋅ F qd (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 aŝ 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][96][97][98] can be simplified greatly 72 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 [K qd (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 3 , 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 Vg(Q) = −Egap + kρ(Q) 2 /2 with Egap = 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 3 , neglecting the residual couplings does not significantly affect the wavepacket [compare panels (a) and (b)], power spectrum I(ω) [panel (c)], or population P ad . 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 3 .
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. of Chemical Physics

B. Induced Renner-Teller effect in hydrogen cyanide
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. 41 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).
Similar to Sec. III A, we simulate the dynamics following an electronic transition from the ground vibrational eigenstate of Vg(ρ, ϕ) = −Egap + E h (ρ) with Egap = 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 P ad , 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 [P ad 1 (t)] ∶= |P ad 1,qd-approx (t)−P ad 1,qd-exact (t)| due to the neglect of the residual couplings is of the same order as the range R P ad 1 ∶= P ad 1,max − P ad 1,min of the population in the whole simulation interval: ϵ res-cpl [P ad 1 (t)]/R P ad 1 = 0.7. Note also that neither P ad 1 (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].

C. Displaced excitation of nitrogen trioxide
Although the excited states of NO 3 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 1 and Q 2 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 residual couplings, however, can be made less important by an improved quasidiabatization. One can reduce the magnitude of the residual couplings from R[F (1) qd (Q)] = 3.8 n.u. to R[F (2) qd (Q)] = 0.5 n.u. by employing the more sophisticated second-order regularized diabatization scheme [45][46][47] obtained by inserting θ (2) (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 P ad 1 (t) due to the neglect of the residual couplings all remain small (see Fig. 6); in particular, quantum fidelity F(t) 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.}

IV. CONCLUSION
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 multistate systems involving multiple conical intersections (including those between three electronic states [104][105][106][107][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.

SUPPLEMENTARY MATERIAL
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).