Numerical Equivalence of Diabatic and Adiabatic Representations in Diatomic Molecules

The (time-independent) Schrödinger equation for atomistic systems is solved by using the adiabatic potential energy curves (PECs) and the associated adiabatic approximation. In cases where interactions between electronic states become important, the associated nonadiabatic effects are taken into account via derivative couplings (DDRs), also known as nonadiabatic couplings (NACs). For diatomic molecules, the corresponding PECs in the adiabatic representation are characterized by avoided crossings. The alternative to the adiabatic approach is the diabatic representation obtained via a unitary transformation of the adiabatic states by minimizing the DDRs. For diatomics, the diabatic representation has zero DDR and nondiagonal diabatic couplings ensue. The two representations are fully equivalent and so should be the rovibronic energies and wave functions, which result from the solution of the corresponding Schrödinger equations. We demonstrate (for the first time) the numerical equivalence between the adiabatic and diabatic rovibronic calculations of diatomic molecules using the ab initio curves of yttrium oxide (YO) and carbon monohydride (CH) as examples of two-state systems, where YO is characterized by a strong NAC, while CH has a strong diabatic coupling. Rovibronic energies and wave functions are computed using a new diabatic module implemented in the variational rovibronic code Duo. We show that it is important to include both the diagonal Born–Oppenheimer correction and nondiagonal DDRs. We also show that the convergence of the vibronic energy calculations can strongly depend on the representation of nuclear motion used and that no one representation is best in all cases.


Introduction
Non-adiabatic effects within the electronic structure of molecules are important for many physical and chemical processes [1][2][3][4][5][6][7] such as when a chemical reaction alter electronic structure, affecting nuclear dynamics.Non-adiabatic processes are also important in astronomy and atmospheric chemistry, where collisions of free radicals and open shell molecules with spatially degenerate electronic states are often seen [8][9][10][11][12] .Modelling electronically non-adiabatic processes has also been effective in explaining the bonding in dications such as BF 2+ 13 and strongly ionic molecules, such as LiF 14 and NaCl 15 , whose 1 Σ + states show non-adiabatic behaviour.
Both the adiabatic and Born-Oppenheimer (BO) approximations assume nuclear dynamics to evolve on single electronic potential energy surfaces (PESs), 8 where no kinetic energy coupling (DDR) to neighbouring electronic states occurs and is generally good for predicting near-equilibrium properties for many molecules 6 .Whilst related, the adiabatic approximation differs from the BO approximation by the addition of the well-known diagonal BO correction (DBOC), introducing mass-dependence into the PECs within the adiabatic representation.The adiabatic approximation then breaks down when electronic states of same symmetry near spatial degeneracy and exhibit an avoided crossing.Neumann and Wigner 16 formalised this as a non-crossing rule for diatomics, showing that potential energy curves (PECs) cannot cross and appear to 'repel' upon approach (see Fig. 1 for example).
Relaxation of the BO and adiabatic approximation is then required to fully encounter the electronically non-adiabatic effects because of the inherent coupling between electronic and nuclear degrees of freedom (DoF) for both the diagonal and non-diagonal terms.
The so called derivative couplings (DDRs) or non-adiabatic couplings (NACs) between states that exhibit avoided crossings arise through the nuclear kinetic energy operator acting on the electronic wavefunctions when the BO approximation is relaxed and corresponds to derivatives in terms of the nuclear coordinate.Computation of DDRs and PESs around the avoided crossing geometry are a major source of computational expense within both quantum chemistry and nuclear motion calculations because of the cusp-like behaviour of the PESs and the singular nature of the DDRs at the geometry of spatial degeneracy 8,[17][18][19] .It is therefore the main focus of many works to explore the property-based diabatisation methods 8,[20][21][22] that transform to a diabatic representation, where DDRs vanish or are reduced and PESs become smooth.For diatomics the smoothness condition of their PECs uniquely defines the unitary transformation to the diabatic representation where NACs (first-order non-diagonal DDR) vanish, PECs are allowed to cross and consequently the molecular properties are smooth, at the cost of introducing off-diagonal diabatic potential couplings (DC).This smoothness is then favourable for nuclear motion calculations since no quantities within the molecular model are singular/cusped making their integration and fitting of analytical forms much simpler.The other method of diabatisation, known as point-diabatisation 14,[22][23][24][25][26][27][28][29][30] , is direct and requires the NAC being obtained ab initio such as through the DDR-procedure 31 where each point can be diabatised without knowledge of the previous one, unlike property-based methods.
Mead and Truhlar 32 showed that a strictly diabatic electronic basis, in which all derivative coupling vanishes, can be defined for a diatomic system.The conditions required to make the first-order NAC vanish are straight forward, however a true diabatic electronic basis only exists when one can remove the second-order (diagonal) derivative coupling simultaneously, which is only possible when considering an isolated state system, allowing one to ignore coupling to other adiabatic states.The adiabatic to diabatic transformation (AtDT) for the N -nuclear-coordinate case up to coupled 4-state systems has been investigated thoroughly by Baer and co-authors since the late 1980s [33][34][35][36][37] .These works develop the so called lineintegral approach in solution to the matrix differential equation that arises when solving for the AtDT which completely reduces the NAC matrix.Their results, albeit from a different angle to this study, are consistent with the results we present.
Despite diabatisation being used routinely to treat the avoided crossings of molecular PESs, there have been very few studies examining the numerical equivalence of adiabatic and diabatic states.This would not only be of value to those who want to benchmark their own nuclear motion codes, but to better understand the roles of each term in the diabatic and adiabatic Hamiltonian.Equivalence refers to the principle that the two representations should yield identical observables, such as energy eigenvalues.
The solution of the nuclear motion Schrödinger equation should not depend on whether the adiabatic or diabatic representations of the electronic states are used 37 .In practice with numerical applications, observables should converge to the same values with increasing accuracy of calculation, e.g by using increasingly larger basis sizes.Equivalency is often assumed but rarely shown.Convergence between the adiabatic and diabatic states have been investigated in only a handful of papers.Zimmerman and George 38 performed numerical convergence tests on adiabatic and diabatic states of the transition probability amplitudes in collisions of collinear atom-diatom systems, where the convergence to equivalence was demonstrated and it was shown that convergence was markedly different with the diabatic representation converging significantly faster.Shi et al. 39 evaluated numerical convergence rates of adiabatic and diabatic energy eigenvalues and eigenfunctions using a sinc-DVR method; equivalency was demonstrated, but this required using a complete adiabatic model and a conical intersection at high energy.The magnitude of the DDR corrections within the adiabatic representation has been studied before such as in the series of papers by Wolniewicz, Dressler and co-workers [40][41][42][43][44][45][46] where excited electronic states of molecular hydrogen and their coupling were studied in detail.The earliest of these studies used the adiabatic approximation but through the series non-adiabatic couplings were introduced and improved for increasing number of excited states and were shown to be essential to produce an accurate spectroscopy (i.e.accurate rovibronic energies and transitions) of the system, confirmed by comparison to experiment.In the later studies the diabatic representation was also shown to provide an accurate description of the nuclear dynamics of H Non-adiabatic interactions are also important for scattering calculations which often assume the equivalence between the adiabatic and diabatic representations 49 .For example, This study aims to show numerical equivalence of the adiabatic and diabatic representation in nuclear motion calculations of rovibronic energies and spectral properties for two selected diatomic systems, represented by two coupled electronic states: yttrium oxide (YO) and carbon mono-hydride (CH) molecules illustrated in Fig. 1.YO shows avoided crossings between the B 2 Σ + , D 2 Σ + and A 2 Π, C 2 Π states as described by Yurchenko et al. 53 .YO has broad scientific interest, having been observed in stellar spectra, [54][55][56][57] found use in solar furnaces 58,59 and magneto-optical traps. 60,61Yttrium oxide is a complex system showing many low-lying electronic states: accurate descriptions of its avoided crossings will be valuable to works in several fields.CH is one of the most studied free radicals 62 because it occurs in such a wide variety of environments: it has been observed in flames 63,64 , solar and 65-67 , stellar spectra 68-70 , spectra of comets 71 , ISM [72][73][74][75] , molecular clouds 76 .

Little and Tennyson
diaba 2 Description of the diabatisation of a two-electronic-

state system
Consider a coupled two-electronic-state system of nuclear (pure vibrational) Schrödinger equations for a diatomic molecule in the adiabatic representation, with the non-adiabatic effects between these two states fully accounted for, as given by (ignoring spin and rotation angular momenta) where r is the distance between the two nuclei and the Born-Huang 2 × 2 Hamiltonian operator is (see, e.g.Varga et al. 30 , Römelt 78 , Yarkony et al. 79 ) Here µ = m 1 m 2 /(m 1 + m 2 ) is the reduced mass, V 12 (r) is the first order DDR, or non-diagonal NAC, given by where ψ a 1 and ψ a 2 are the adiabatic electronic wavefunctions, and K(r) is the diagonal DDR component given by Furthermore, ℏ 2 2µ K(r) is the well-known DBOC 80 .The derivative coupling K(r) is related to the second DDR W In conjunction with Eqs.(4-6) and the results by Baer 37 , N. Mabrouk and Berriche 83 , Smith 84 , a simple and powerful expression for the matrix element of the diagonal DDR term K for the coupled two-electronic state problem is obtained: = dW (1) 12 dr − W (2) .
A diabatic representation of a two-state system can be introduced via a unitary transformation U(r) of the adiabatic electronic wavefunction vector ψ a = (ψ a 1 , ψ a 2 ) T in which the 1 st order DDR vanishes and PECs and other molecular properties become smooth at the cost of introducing an off-diagonal potential energy coupling, termed a diabatic coupling (DC), between the non-adiabatically interacting electronic states 17,18,85 .The unitary 2 × 2 matrix U(r) is given by where the mixing angle β(r) is obtained by integrating NAC as follows 8,86-88 where r 0 is a reference geometry and is usually chosen as such that one can define a physical condition which ensures the mixing angle to equal π/4 at the crossing point r c .It can be also shown that for the diatomic one-dimensional case the transformation to a strict diabatic basis is unique and that W (1) 12 vanishes upon the diabatisation together with K(r) (see Eq. ( 7)).Similar to the work by Köppel et al. 89 who developed a Hamiltonian for the two-coupled electronic state problem, we develop theory for the diabatic and adiabatic electronic potential energy curves for the coupled two-electronic states in question.The corresponding two-electronic-state Born-Huang Hamiltonian operator Ĥd then becomes where the diabatic potential energy functions V d 1 (r) and V d 2 (r) and the DC function are given by The goal of this work is to demonstrate the equivalency of the adiabatic and diabatic representations when solving the nuclear motion diatomic (eigenvalue) problem.To this end we aim to construct, solve and compare the eigensolutions of model diatomic systems in the adiabatic and diabatic representations.
If the adiabatic representation of an isolated two-electronic state diatomic is fully defined by the three functions V a 1 (r), V a 2 (r) and W 12 (r) in Eq. ( 1), in turn, the diabatic representation is fully defined by the three functions V d 1 (r), V d 2 (r) and V d 12 (r) in Eq (10).In fact, both transformations can be fully described by a combination of any three functions from 12 (r), V d 1 (r), V d 2 (r) and V d 12 (r).For this study we choose V d 1 (r), 12 (r).The diabatic PECs V d 1 (r), V d 2 (r) are expected to have smooth shapes by construction and are easy to parameterize, which explains our choice, while W 12 (r) has also a rather simple, easy-to-parameterize cusp-like shape 8,14,[17][18][19] as will be shown below.
The other three functions are constructed from V d 1 (r), V d 2 (r) and W (1) 12 (r) as follows.
We first define β(r) via Eq.(9).By applying the inverse transformation U † to the potential matrix V d (r) in Eq. (11), we arrive at the following condition for the off-diagonal element of the adiabatic potential matrix which is required to be zero since V a (r) = UV d (r)U † in Eq. ( 1) is diagonal by definition.
Hence, we can rearrange it for the DC to get The adiabatic functions V a 1 (r) are V a 2 (r) can be then constructed as eigenvalues of the diabatic potential energy matrix (second term in Eq. ( 10)): or, equivalently, via the inverse unitary transformation U:

Spectroscopic Models
As an illustration, two model two-state electronic systems are used, YO and CH with their diabatic and adiabatic curves shown in Fig. 1 and introduced in detail in the following.

YO spectroscopic model
As an example of a two-state system with narrow, coupled bound electronic curves, we choose the ab initio PEC curves of the B 2 Σ + and D 2 Σ + states of yttrium oxide from Smirnov et al. 90 with the NAC from Yurchenko et al. 53 .
We use a Morse oscillator function as a simple model for the diabatic B 2 Σ + and D 2 Σ + PECs of YO as given by where A e is a dissociation asymptote, A e − V (r e ) is the dissociation energy and r e is an equilibrium distance of the PEC.The NAC of YO can be efficiently described by a Lorentzian function: where γ is the corresponding half-width-at-half-maximum (HWHM), while r c is its center, corresponding to the crossing point of diabatic curves.These PECs and NAC are illustrated in Fig. 2. The parameters defining these curves are listed in Table 1, which were obtained by fitting to the corresponding ab initio data.
For the Lorenztian as a NAC, Eq. ( 9) is easily integrable to give the transformation angle β(r): where r c is obtained as the crossing point between the PECs, and the ± sign refers to the path integral when r < r c and r c < r respectively.
The adiabatic curves obtained using Eqs.(14,15) and the DC curve obtained using Eq. ( 13) are shown in Fig. 2. The value of the crossing point r c is obtained as numeri- and is listed in Table 1.The derivative coupling K in the diagonal matrix element of the adiabatic kinetic energy operator in Eq. ( 1) is simply defined by 2 according to Eq. ( 7).All the corresponding curves are programmed in Duo analytically and are provided on a grid of 1000 equidistant bond lengths as part of the supplementary material.
Table 1: The molecular parameters defining the YO spectroscopic model T e , cm −1 20700.020400.0

CH spectroscopic model
The spectroscopic model for CH, with curves illustrated in Fig. 1 (right panel), is constructed to mimic the ab initio curves of C 1 Σ + and 2 1 Σ + by van Dishoeck 91 .The C 1 Σ + state has a bound shape with a well of about 16700 cm −1 (2.0705 eV), which we model using a Morse oscillator function in Eq. ( 17).The 2 1 Σ + state is repulsive, with the dissociation energy lower than that of C 1 Σ + by about 10000 cm −1 .We chose to model the 2 1 Σ + PEC using the following form: The corresponding NAC between C 1 Σ + and 2 1 Σ + of CH from van Dishoeck 91 is modelled using a two-parameter Lorentzian function in Eq. ( 18).All parameters defining the CH spectroscopic model are listed in Table 2.As above, the value of the crossing point r c is obtained as a numerical solution of Table 2: The molecular parameters defining the CH diabatic spectroscopic model T e , cm −1 32500.00.0 r e , Å 1.12 1.12 b, Å−1 2.5 1.968A e , cm −1 49200.029374.039220.0  1) or ( 10) are extended with the rotation-spin-electronic contribution as follows (see Yurchenko et al. 77 for details of the approach used): where the rotational angular momentum operator R is replaced with Here Ĵ, Ŝ, L are the total, spin and electronic angular momenta, respectively.We then solve the aforementioned rovibronic Schrödinger systems for YO and CH variationally on the Hund's case (a) basis using the Duo program, 77 which has been extended as part of this work to include the adiabatic and diabatic effects.The spectroscopic models of CH and YO are provided in the form of the Duo input files both in the diabatic and adiabatic representations as the supplementary material.
Duo uses the numerical sinc-DVR method 92,93 to solve the Schrödinger systems for the curves defined either on a grid or as analytic functions.For the analytic representations as above, the corresponding functions are mapped on a grid of sinc-DVR points.For a grid input, cubic splines are used.The Duo kinetic energy has been extended to include the first derivative component required for implementation of the NAC, also using the sinc-DVR representation 94 .The DBOC terms can be either provided as input or generated from the NAC using Eq. ( 3).In order to facilitate numerically exact equivalency of the diabatic and adiabatic representations Duo calculations, Eqs.(13, 14 and 15) are provided are used for constricting V d 12 , V a 1 (r) and V a 2 (r), respectively, from V d 1 (r), V a 2 (r) and β(r).

The YO solution
We first find the vibronic (J = 0.5) energies of the coupled B 2 Σ + and D 2 Σ + systems in the adiabatic and diabatic representations as accurately as possible in order to establish a baseline and also to demonstrate the equivalency of the two representations.Even though we know that the diabatic and adiabatic solutions should be equivalent (i.e.identical within the calculation error), this is always a subject to the convergence or other numerical limitations.
For example, Duo uses a PEC-adapted vibrational basis set constructed by solving the pure vibrational problem, which will be different depending on the representation, diabatic or adiabatic, and thus will influence the convergence.The corresponding YO model curves are shown in Fig. 2, where the DBOC coupling K is included into the adiabatic PECs for clarity.
There is a striking difference between the two models, with a large spike in the middle of the adiabatic PECs, yet we expect them to give the same eigenvalues and eigenfunctions.
A selected set of rovibronic energy term values (J = 0.5) computed using the two methods is listed in in the adiabatic model, the diagonal DDR is switched off (K = 0), but the NAC is kept in; (A3) in the diabatic model, the diabatic coupling is set to zero (V 12 = 0).The effects of these approximations on the calculated energies of YO (J = 0.5) are also shown in Table 3.
For the adiabatic model, the omission of K (A2) has the overall largest impact especially on the B 2 Σ + term values.The omission of V 12 from the diabatic model (A3) appears to be less damaging than other two approximations.It is clear, however, that any degradation of theory leads to large errors, unacceptable for high-resolution applications.This is in fact the main conclusion of this work that the impact of dropping any of non-adiabatic corrections from the model describing a system with crossing has to be always investigated.
Out of the two representations, the adiabatic model is usually considered to be more complex to work with.Its curves have complex shapes with the model being very sensitive to the mutual consistency of the curves V a 1 , V a 2 and W 12 around the crossing point.The disadvantage of the diabatic representation is that it does not come out as a solution of the (adiabatic) electronic structure calculations directly and needs to be constructed either through a diabatisation approach 8,14,20-30 or approximated.

Eigenfunctions and Reduced Density
It is instructive to compare the eigenfunctions φ J,τ i (r) of the adiabatic and diabatic solutions and different approximations.To this end we form reduced radial densities of the eigen-state in question.The eigenfunctions φ J,τ i utilised by Duo are expanded in the basis set |n⟩, where N is the basis size and C J,τ i,n are the expansion coefficients used to assign quantum numbers by largest contributions.|n⟩ denotes the full basis: |n⟩ = |st, J, Ω, Λ, S, Σ, v⟩ where 'st' is the electronic state, S is the electron spin angular momentum, v is the vibrational  quantum number, and Λ, Σ and Ω are the projections of electron orbital, spin, and total angular momentum along the internuclear axis respectively.The reduced radial density ρ J,τ i (r) is then given by where |k⟩ = |st, J, Ω, Λ, S, Σ⟩ and χ v (r) are the vibrational wavefunctions.The reduced density states are probability density functions over bond length averaged over all quantum numbers in |n⟩.This is an efficient way of examining the behaviour of the wavefunctions without looking in a hyperdimensional space defined by quantum numbers |n⟩.
Figure 3 shows selected reduced radial state densities of YO computed using different representations and approximations.As expected from our energy comparisons, the diabatic and adiabatic representations produce identical results, whereas the reduced densities quickly deviate when omitting the NAC and/or K corrections.Again, it appears that the adiabatic representation with approximations is almost better when the DDRs are completely omitted rather than omitting only one, at least concerning the lower energy levels.

Adiabatic and diabatic solutions for CH
We now turn to a slightly different system of the C 1 Σ + and 2 1 Σ + states of 12 CH shown in Fig. 4. Adiabatically, these states have a large separation and a broad NAC.In contrast to YO, there is no spike-type contribution from the DBOC-term K to the adiabatic PECs of CH.Diabatically, the system consists of a bound and a repulsive state with a crossing at large distance and high energy which therefore should not influence the lower rovibronic states of C 2 Π significantly.Regardless of the representation used, the region above the first dissociation channel (39220.0cm −1 ) is heavily (pre-)dissociated and should contain both (pre-)dissociative and continuum states.Duo is capable of finding both the bound and continuum eigensolutions.While the bound wavefunctions satisfy the standard boundary condition to decay at large and short distances, the continuum wavefunctions can be also computed with the sinc-DVR method used by Duo and satisfy the boundary condition of vanishing exactly at the simulation box borders (together with its first derivatives), see Pezzella et al. 95 .For the analysis we separate the (quasi-)bound and the continuum states by checking the character of the wavefunctions at the 'right' border r max : while the continuum states tend to oscillate at r → ∞ with a non-zero density around r max 96 where the bound state vanish completely.
The resulting energy term values of the bound states are listed in Table 4 for all five cases, including non-adiabatic and diabatic couplings are considered as in the YO example.
The full diabatic and adiabatic (bound) C 1 Σ + energies are fully equivalent within 10 −6 cm −1 (here shown up to the second decimal point).However any degradation of the theory leads to drastic changes in the topology of the system and hence in the calculated rovibronic energies of the C 1 Σ + state, with the error quickly deteriorating already for v = 2.For example, removing the DC term, the diabatic solution becomes meaningless with lots of non-physical bound states above the first dissociation channel, non-existent in the case of the full treatment.A similar effect is caused by the omission of the derivative couplings from the adiabatic pictures with bound spurious 2 1 Σ + states produced by the adiabatically bound PEC 2 1 Σ + (see Fig. 4).Although the omission of the K(r) term from the adiabatic solution seems harmless for the topology of the corresponding PECs, even this case leads to a spurious vibrational 2 1 Σ + (v = 0) state.Therefore the conclusion is that every non-adiabatic term should be considered important, unless proven otherwise.
The corresponding reduced densities for some lower lying bound states of CH (C 1 Σ + , J = 0.5) are shown in Fig. 5 (n = 1, 2, 3).We see that the low lying vibronic states of C 2 Π are largely unaffected by the omission of the DDRs or DC since they are energetically well separated from the region of non-adiabatic interaction, in this case occuring near dissociation.However, the reduced densities of the 2 1 Σ + state (n = 4) quickly diverge when removing the NAC and/or K correction.The 2 1 Σ + is adiabatically bound and diabatically unbound, where this drastic difference is seen in the reduced densities of Fig. 5, and correspond to energy levels which arise from PECs of very different character.For example, in the diabatic case where the DC is omitted, the n = 4 state corresponds to the bound C 1 Σ + (J = 0.5, +, v = 3) state whereas in the adiabatic A1 and A2 cases the n = 4 bound state corresponds to the bound 2 1 Σ + (0.5, +, v = 0) state.In the cases where the DDRs and DC are fully accounted for no fourth bound state exists since the couplings will push it into the quasi-bound region about the adiabatic potential hump of the C 1 Σ + state.This quasi-bound nature begins to show itself in the reduced density of the adiabatic case with K = 0 where small oscillations propagating to the right simulation border at 4 Å are seen.

Continuum solution of CH: photo-absorption spectra
In order to illustrate the equivalence of the continuum solution involving in the repulsive 2 1 Σ + of CH, we model a photo-absorption spectrum X 1 Π → C 1 Σ + /2 1 Σ + , where we follow the recipe from Pezzella et al. 97 and Tennyson et al. 98 .For the X 1 Π state, we use the same Morse function representation in Eq. ( 17) with the parameters listed in Table 2.For the transition electric dipole moments μX,C = ⟨X 1 Π| µ |C 1 Σ + ⟩ and μX,2 = ⟨X 1 Π| µ |2 1 Σ + ⟩ of CH we adopt the ab initio curves by van Dishoeck 91 with an approximate model using the following function: where ξ p is the Šurkus 99 variable given by: The parameters defining the diabatic transition dipole moment (TDM) functions are listed in Table 5.The adiabatic TDM curves are obtained through the unitary transformation U (r): Table 4: The rovibronic (J = 0.5, 1.5 and 2.5) bound energy term values (cm −1 ) of the C 1 Σ + (C) and 2 1 Σ + (2) systems of CH computed using the adiabatic and diabatic representations.
Table 5: The molecular parameters defining the CH diabatic transition dipole moment functions  using the continuum solution of the coupled C 1 Σ + /2 1 Σ + system from the bound states of X 1 Π, for the diabatic and adiabatic models.We used the box of 60 Å and 1600 sinc-DVR points.For the cross sections, a Gaussian line profile of the half-width-at-half-maximum (HWHM) of 50 cm −1 was used to redistribute the absorption intensities between the 'discrete'

Convergence
Since Duo uses a solution of the J = 0 uncoupled vibrational problem to form its vibrational basis set functions ψ v (r), and these model problems are hugely different depending on the representation, one can also expect that the convergence of the eigensolution to be impacted by the choice of the representation.
Here we test the convergence of the J = 0.5 energy levels of our simplified YO and CH models in the diabatic and adiabatic representations where all non-adiabatic effects encountered.Figure 9  Tests comparing the convergence rates for vibrational energies of higher J resulted in the same conclusions as above for the J = 0.5 case.
This shows that there is not one representation that rules over the other, it depends on the character of the avoided crossing, specifically in its position, the shape of the potentials approaching the crossing, and the separation of the adiabatic PECs.It is therefore important to consider the system of study before choosing a representation, where all corrections must be included.is plotted as a function of vibrational basis size.A constant grid size of G= 3001, 4001 points for the sinc-DVR basis set was used for the YO and CH states, respectively.We see that the diabatically computed energies for YO converge much faster than the adiabatic ones, whereas for CH the opposite is true.

Conclusion
A demonstration of the equivalency of the diabatic and adiabatic representations for two model diatomic systems, bound electronic B 2 Σ + and D 2 Σ + states of YO and abound/repulsive electronic system C 1 Σ + and 2 1 Σ + of CH, is presented.Both representations should be equivalent by construction, but we explicitly show this within nuclear motion calculations through comparison of the rovibronic energies and wavefunctions.The importance of different non-adiabatic couplings in the molecular Hamiltonian are investigated, such as how the rovibronic energies and wavefunctions change when the NAC, the second derivative coupling or the diabatic coupling vanish.
We present a transformation from the adiabatic to the strict diabatic basis for an isolated two-electronic state diatomic system.Each representation is defined by three functions, the adiabatic representation is given by two avoiding PECs and their corresponding NAC whereas the diabatic picture is analogously defined by two diabatic PECs and a DC, all of which are related to each other through the mixing angle.Because of this, any three of the aforementioned quantities can be used to fully reconstruct either the adiabatic or diabatic representation.We demonstrate that the choice of two diabatic PECs and a NAC provides an easily paramaterisable and powerful way to define the two-level problem.In the case of the diabatic PECs, they can be modelled easily by Morse oscillators, and the NAC is easily modelled using a Lorentzian.
We show that omission of any of the non-adiabatic terms lead to significant changes in the spectral properties of these systems, unsatisfactory especially for high-resolution applications.Even the diagonal derivative coupling, often omitted in practical applications, is shown to be of central importance when achieving equivalency.
We also show that the choice of a preferable representation, diabatic or adiabatic, is not the same for all systems.For cases where the NAC is small (large DC) then the adiabatic representation shows fast convergence of rovibronic energy levels.However, for cases where the NAC is large (small DC) then the diabatic representation converges rovibronic energies with very small basis sets where large ones are required for the corresponding adiabatic representation.
We used simplified approximated functions to model different diabatic and adiabatic curves for the purpose of facilitating the comparison and demonstration of the equivalency, as well as to simplify the debugging process.In fact, our program Duo uses numerically defined curves either provided as grids of r-dependent values or generated from analytic input functions as used here.For the convenience of the reader all curves from this work are provided both in the analytic and numerical representations as ASCII files, which are also Duo input files.As we demonstrated, the models provide exact equivalency of the diabatic and adiabatic salutations and therefore can be used as benchmark for similar programs.At It would be interesting to develop and apply a similar methodology for polyatomic molecules, where the derivative couplings cannot be fully transformed away.The exact equivalence of the two representations should be still possible to demonstrate numerically even for a quasi-diabatic transformation.This work for triatomic molecules is currently underway.In the present diatomic study we show that exclusion of the DDR couplings can lead to differences on the order of magnitude of 10s-100s of cm −1 in the energy wavenumbers, reinforcing the need for a careful error budget of all the approximations made when using them in high-resolution spectroscopic applications.

Associated Content
All the DDR, potential energy, and DC curves are programmed in Duo analytically and are provided on a grid of 1000 equidistant bond lengths as part of the supplementary material.
The spectroscopic models of CH and YO are also provided in the form of Duo input files both in the diabatic and adiabatic representations as part of the supplementary material.

Conditions for a Strictly Diabatic Basis: Nuclear Motion Part
Let us begin with the strict diabatic basis, a frame where both the off-diagonal DDR and DBOC couplings vanish simultaneously, which was shown to be possible for the diatom by Mead and Truhlar 32 .Within this diabatic representation the Hamiltonian comprises of a purely diagonal kinetic energy operator and an electronic potential matrix with non-zero off-diagonal coupling elements analogous to the non-adiabatic couplings (DDRs) within the adiabatic representation, called diabatic couplings 17,18 (DCs).All PECs are allowed to cross in this frame and property operator curves (e.g.SOCs) are smooth.We now wish to study the diabatic and adiabatic Hamiltonians in a similar vein to the work by Köppel et al. 89 who develop a Hamiltonian for the two-coupled electronic state problem.With this, the diabatic coupled-two-electronic state Born-Huang Hamiltonian takes the following matrix form where V d 1,2 is the DC (potential) coupling the two electronic states in question, |φ d i ⟩ are the diabatic nuclear basis functions, r is the diatom nuclear bond coordinate, and the superscripts denote the adiabatic 'a' or diabatic 'd' representation.To find the conditions required to achieve this strict diabatic basis, we transform back to the adiabatic representation via a inverse unitary transformation U † where U transforms the molecular system from the adiabatic to diabatic (AtDT) representation -the well known ATD transformation 17,[33][34][35][36][37]79,100 -via the following rotation We define our transformation to the diabatic representation through the transform of the coupled PECs (see text) Hence the diabatising transformation can be thought to either mix operator matrix elements or mix the adiabatic basis functions through a rotation of angle β(r This rotation is called the mixing angle 8,79,[86][87][88]100 and is a function of the diatom bond length. We cn now transform the diabatic Hamiltonian to the adiabatic representation through action of U on the nuclear kinetic energy operator and then enclose with adiabatic nuclear basis to find the adiabatic matrix elements T (a) α,β , we have in atomic units, Before evaluating the above matrix elements, we introduce a minus sign and invert the direction of one of the first derivative components of the second derivative via a Hermitian conjugation which can be done under integration of nuclear coordinates.We show this explicitly via the following integration by parts in wavefunction notation, Since the wavefunction vanishes for |r| → ∞ then we have φ a * α U d dr U † φ a β ∞ −∞ = 0. Thus, we are left with We will see that this inversion of the direction of derivatives will ultimately lead to the adiabatic Hamiltonian being in its Hermitian form and will be convenient for us to find the conditions of a strictly diabatic basis.Evaluating the derivatives using the product rule and expanding the equation we find the following matrix elements, T since U U † is the identity.This infers that the products dU dr U † and U dU † dr are equal to a skew-symmetric matrix since S + S † = 0. Hence, We can now see that the derivative matrix product dU dr dU † dr is given in terms of the skewsymmetric matrix via For 2 × 2 unitary matrices, the corresponding skew symmetric matrix S only has one unique element, ζ, and the squared matrix S 2 has every diagonal element equal to −ζ 2 .Thus, for a 2-electronic state problem, the unitary matrix product dU dr dU † dr is diagonal.Through a similar argument, the matrix product U dU † dr is a skew-symmetric matrix, with zero in the diagonal.We can thence write the kinetic energy matrix element T (a) α,β as The first term contains the diagonal matrix elements and the second term are the off-diagonal terms.As we will see in the next section, U dU † dr is equivalent to the first-order NAC matrix, W (1) , and defines the condition to make this NAC vanish upon action of the AtDT.The first-order matrix differential equation U dU † dr = W (1) has been investigated thoroughly by Baer and co-authors since the late 1980s [33][34][35][36][37] .These works develop the so called line-integral approach in solution to the above matrix differential equation, which for a two-level system, is easily achieved by noticing the common solution to the first order ODE U = e − W (1) dr .
It turns out this solution is exact for the diatomic two-coupled electronic state problem and arises because W (1) and its integral commute.One can then manipulate the taylor series expansion for U = e − W (1) dr to show that it is indeed equal to a rotation matrix of the mixing angle β.To this end, let us now insert the rotation matrix form for U into the above expression for the kinetic energy matrix elements, noting the following matrix multiplications Since at this point the electronic DoF have been already integrated over, we can identify the scalar terms in the above Hamiltonian with the DDRs one would get from ab initio calculations.The only scalar term that arises is the derivative of the mixing angle, and therefore the condition to transform to the strict diabatic basis depends solely on the transformation through dβ dr , the so called Abelian case as discussed by Baer 37 .We know that the DDRs will consist of a first DDR, a.k.a. the NAC, which we call generally W (1) here which is antihermitian, and a second DDR diagonal coupling element which we call K. Thus, the condidtions for a strictly diabatic basis can be summarised by which we will show in the next section to be a sensible set of condiditions, thus providing a way to form the strict diabatic basis for diatomic molecules in the isolated two-electronic state problem.

The Non-Adiabatic Derivative Couplings: Electronic Structure Part
We now turn to the electronic structure part of our discussion on adiabatic and diabatic representations for the diatomic molecular Hamiltonian.Let us begin with the electronic matrix elements of the nuclear kinetic energy operator, or second derivative with respect to the nuclear bond length where |ψ a ⟩ are the adiabatic electronic basis wavefunctions and W α,β are the second DDR derivative coupling elements.By evaluating these matrix elements we are relaxing the Born-Oppenheimer approximation which would assume these elements to be zero.The NAC, or 1 st DDR, is well known and defined as the matrix elements of the first derivative where the NAC is anti-Hermitian and can be investigated using the Hellman-Fynman theorum.To this end, consider the following eigenvalue equation Where in the last part we used the Hermiticity of the Hamiltonian to act Ĥ(a) on the bra yielding us E α .We now have, Hence, if states α and β are sufficiently well seperated, i.e. |E β − E α | ≫ 1, then the DDR matrix elements are small W α,β ≪ 1.Multiple studies 37,81,82 have used the following identity to relate W (2) and W (1) via where the DBOC K α,β = dψ a α dr dψ a β dr ̸ = 0 for α = β.Historically, the use of the above identity was to avoid the cumbersome computation of W (2) by only having to evaluate W (1) and its derivative, and is identical to the g-, h-, and k-notation by Lengsfield and Yarkony 81 , Saxe and Yarkony 82 with h = dg dr − k.The above relation shows the second DDR to have both diagonal and off-diagonal components and can be compared to the scalar terms in Eq. (29).
We can simplify further our expressions for the DDRs by introducing a resolution of the identity between the bra and ket of the DBOC K ρ,ρ within the adiabatic basis, yielding where the summation is over all adiabatic states and in the last line we inserted the Hellman-Feynman relation in Eq. (30).We see that for the coupled two-electronic state system, states |ψ a 1 ⟩ and |ψ a 2 ⟩, the NAC elements coupling other states will be reduced by a factor of 1 (Eρ−Eκ)(Eκ−Eρ) , for sufficiently well separated states from the coupled system |ψ a 1 ⟩ and |ψ a 2 ⟩ we can truncate the summation to over these two states only.We now have which is diagonal and equal to minus the NAC squared.We now see that all DDRs are completely defined in terms of the 1 st DDR, or NAC, only.Knowing we will want to find matrix elements within the nuclear basis later, we rewrite Eq.( 31) in matrix form The condition that the 2 nd DDR must equal the square of the NAC in order to completely remove the DDR couplings is then sensible for a coupled two-electronic state system which is energetically well-separated from other adiabatic states, as shown in Eq.( 32) under Hellman-Fynman formalism.This of course only works when the two-level system is effectively isolated, and breaks down upon approach of closely lying adiabatic states with common symmetries. (

Figure 1 :
Figure 1: Illustration of the [D 2 Σ + , B 2 Σ + ] and [C 2 Σ + , 2 2 Σ + ] avoided crossing systems ([black, blue] lines) for the YO and CH diatomics, respectively, which we use to perform tests on the adiabatic and diabatic equivalence.The top panels show the adiabats (solid lines) and diabats (dashed lines).The bottom panels show the corresponding NAC and DC (in the units of Å and cm −1 ) of the transformations.
Solving the rovibronic Schrödinger equations for CH and YO Both the CH and YO doublet systems represent open-shell molecules.Towards a complete rovibronic solution, the pure vibrational Hamiltonian operator in Eqs. (

Figure 2 :
Figure 2: Full adiabatic (left) and diabatic (right) models of the B 2 Σ + and D 2 Σ + systems of YO.The top panels show the PECs, where the adiabatic PECs include the diagonal DDR correction αK and α = h/(8π 2 cµ).The bottom panels show the corresponding coupling curves, NAC (left) and DC (right).

Figure 3 :
Figure 3: YO reduced density states for the lowest 5 bound levels with n being the energy enumerator given in Table3.These reduced densities are illustrated and computed using different levels of theory: diabatic representation with DC (blue dotted); diabatic model with the DC turned off (magenta, A3); adiabatic representation with both the NAC and K correction included (lime green); adiabatic representation with NAC only (orange, A2); adiabatic representation with no correction (red, A1).

Figure 4 :
Figure 4: Full adiabatic (left) and diabatic (right) models of the C 1 Σ + and 2 1 Σ + systems of CH.The top panels show the PECs, where the adiabatic PECs include the diagonal DDR correction αK, where α = h/(8π 2 cµ).The bottom panels show the corresponding coupling curves, NAC (left) and DC (right).

Figure 5 :
Figure 5: CH reduced density states for the lowest four bound rovibronic levels with n being the energy enumerator given by the row number in Table4.Different levels of theory are used to compute these reduced densities and are illustrated: diabatic representation with DC (blue dotted); diabatic model with the DC turned off (magenta, A3); adiabatic representation with both the NAC and K(r) correction included (lime green); adiabatic representation with NAC only (orange, A2); adiabatic representation with no correction (red, A1).

Figure 7
Figure7shows a photo-absorption spectrum of CH at T = 300 K computed with Duo

Figure 6 :
Figure 6: Adiabatic (left) and diabatic (right) models of the photo-absorption system of X 1 Π → C 1 Σ + /2 1 Σ + of CH.The top panels show the PECs, adiabatic and diabatic, while the bottom panels show the corresponding transition dipole moment curves.

Figure 7 :
Figure 7: Photo-absorption spectra of CH at T = 300 K.The no-approximation case is shown with the blue line; the NAC=0 case is shown with the red line and the black line shows spectrum with all DDRS set to zero.

Figure 8 :
Figure 8: Reduced density of the continuum state corresponding to an energy of hc • 38183.6576cm −1 .Its transition with the X 1 Π(J = 1.5, f, v = 0) state is positioned at the peak in the spectra of Fig.(7).The reduced density state is illustrated and computed using different levels of theory: diabatic representation with DC (blue dotted); diabatic model with the DC turned off (magenta, A3); adiabatic representation with both the NAC and K correction included (lime green); adiabatic representation with NAC only (orange, A2); adiabatic representation with no correction (red, A1).
illustrates the convergence of the lowest 20 J = 0.5 energies of YO and the n = 5 state of CH (C 2 Σ + (J = 0.5, ±)) where the difference of the i-th level Ẽi to its converged value Ẽcvg i is plotted as a function of vibrational basis size.The two systems show contrasting results.The diabatically computed YO (D 2 Σ + ) energies converge very quickly for basis sizes of ∼ 25 whereas within the adiabatic representation a much larger basis set of ∼ 250 was required to achieve convergence.For CH (C 2 Σ + ) the adiabatic energies initially converge faster but the diabatic energies eventually converge to within 10 −6 cm −1 for a basis size of ∼ 25 as opposed to ∼ 42 for the adiabatic energies.

Figure 9 :
Figure 9: Convergence of the lowest 20 vibrational J = 0 energies of the D 2 Σ + state of YO (left) and the C 2 Σ + (v = 0, e/f ) state of CH (right) are plotted, where the difference of the i-th vibrational level E i to its converged value E cvg i the same time, Duo provides an efficient platform to test different aspects of diabatisations in diatomic calculations, including testing different approximations.Duo is an open-access, with an extended online manual and many examples.

2 ,|I d 2 dr 2 Ĥ
where σ = iσ and σ is the second Pauli matrix.We can now simplify our expression for the adiabatic kinetic energy matrix elements by use of the above relations, giving us Noting that the adiabatising transformation diagonalises the electronic Hamiltonian by construction, performing a Hermitian conjugation of one of the first derivative components in the second term above under integration of the nuclear coordinated to write dφ a α |φ a β ⟩, we can write the adiabatic 2×2 Born-Huang Hamiltonian as, going back to SI units by introduction of the ℏ 2 /2µ kinetic energy factor, is not clear if it is Hermitian (see next section) and off-diagonal terms to the nuclear kinetic energy operator are introduced, which will correspond to the non-adiabatic derivative couplings, or DDRs, in this representation.Here the direction of the derivatives are still shown, and is how we implement the matrix elements within our nuclear motion code Duo.We can again introduce a minus sign and redirect the derivative through Hermitian conjugate to arrive at the following

W
In preperation to finding the matrix elements of the kinetic energy operator in nuclear basis, one would left and right multiply T by the Born-Huang wavefunction |Ψ a tot ⟩ = |ψ a ⟩ |φ a ⟩.However, let us consider only the electronic part and use the relations we have derived in this section, we have⟨ψ a α | T (a) |ψ a β ⟩ = −Where we have explicitly shown the direction of the derivative.Knowing we will take matrix elements within the nuclear basis we now push through the first derivative on the first-order the integration of nuclear DoF we can also perform a Hermitian conjugation on the − → .With this, we now have ⟨ψ a α | T (a) |ψ a β ⟩ = − and we used the antihermiticity of W (1) α,β to introduce another minus sign on the element ⟨ψ a 2 | T (a) |ψ a 1 ⟩.Comparison with the adiabatic molecular Hamiltonian inEq.(28)yields the final conditions for transforming to a strictly diabatic basis to be W 2 , but comparisons between the adiabatic and diabatic representations were not shown.Additionally, DDR corrections were studied with respect to the computed rovibrational energies of H + 2 and D + 2 by Jaquet and Kutzelnigg 47 and later by Jaquet 48 on the H + 2 , H + 3 , and H 2 systems.It is therefore expected that DDR contributions are critical for the accurate determination of the energies of small hydrogen bearing molecules.

Table 3 .
The energies are indeed identical (within 2.5 × 10 −5 cm −1 ), but the approximate quantum state labels as assigned by Duo are very different.Duo assigns quantum labels via the largest contribution from the corresponding basis sets, which in both cases are very different and so are their state interpretations, in which case we compare states of matching energy enumerator n.
Having established the numerical equivalence, we can now investigate the importance of different non-adiabatic couplings for the YO model.Three approximations are considered here: (A1) in the adiabatic model, both DDR terms are switched off (W

Table 3 :
The rovibronic (J = 0.5) energy term values (cm −1 ) of the B 2 Σ + (B) and D 2 Σ + (D) systems of YO computed using the adiabatic and diabatic representations.The energies are listed relative to the lowest J = 0.5 state.