Born-Oppenheimer potentials for $\Pi$, $\Delta$, and $\Phi$ states of the hydrogen molecule

We report on accurate variational calculations of the Born-Oppenheimer potential for excited states of the hydrogen molecule with $\Pi$, $\Delta$, and $\Phi$ symmetries. The obtained potential energy curves reach the relative precision of $10^{-9}$ or better along internuclear distances of 0.01 -- 20 au. Calculations rely on recursive evaluation of two-center two-electron molecular integrals with exponential functions in arbitrary precision arithmetics. Our results for most of the states are the first ever reported, and for the previously calculated states constitute an improvement by several orders of magnitude.


Introduction
The hydrogen molecule is one of the most intensively studied two-electron system. Its rotational and vibrational energy levels in the ground electronic state are known with 10 −9 relative accuracy, by including nonadiabatic effects and relativistic and quantum electrodynamic corrections up to the order of α 6 m, where α ∼ 1/137 is the finestructure constant. Theoretical predictions agree very well the recent high-precision measurements [1][2][3], including those of the dissociation energy, which reach the accuracy of 10 −4 cm −1 [4][5][6]. Similar if not better experimental accuracy is achieved for the transitions to the excited electronic states, but theoretical predictions are far less accurate here.
In our recent work we have performed calculations of Born-Oppenheimer (BO) potentials for excited states of all n Σ + symmetries with n ≤ 7, for which we obtained the accuracy exceeding all the previously known results by 3-4 orders of magnitude. In this work we extend those calculations to Π, ∆, and Φ excited states with the relative accuracy of at least 10 −9 for internuclear distances of 0.01 -20 au. Some of these states have already been investigated in the literature. Namely, the most significant and accurate results at that time were obtained primarily in a series of works by Ko los, Wolniewicz, and Rychlewski [7][8][9][10][11], using the explicitly correlated exponential basis, called Ko los-Wolniewicz (KW) basis functions [12]. Their results, however, disagreed with experimental values for ∆ states. Jeziorski et al. pointed out in Ref. [13] that michal.silkowski@fuw.edu.pl the angular factor for these states in Ref. [11] was incomplete. This was verified and later corrected by Wolniewicz in Ref. [14]. It lowered BO curves by a few cm −1 for singlet and less than 1 cm −1 for triplet ∆ g states and resolved discrepancies with the experimental data available at those times.
Another method, based on explicitly correlated gaussian (ECG) functions, was employed in Ref. [15,16] to obtain a potential curve for C 1 Π u (n = 1). The most recent contributions to the ab-initio Potential Energy Curves (PEC) of H 2 consist of a Free-Complement local Schrödinger equation method developed by Nakatsuji and collaborators [17][18][19][20][21] with potential curves for Σ, Π, ∆, and Φ states. In this work we not only improve all these previous results by a few orders of magnitude, but also calculate many higher excited states of Π, ∆, and Φ symmetries and finally draw attention to possible further applications of Kolos-Wolniewicz basis.

Method
Our variational calculations utilize explicitly correlated exponential functions with polynomial dependence on interparticle distances of the form [12], where η i and ξ i are proportional to confocal elliptic coordinates and are given by η i = r iA − r iB , ξ i = r iA + r iB with i enumerating electrons and real y, x, u, w nonlinear parameters subject to variational minimization. By {n} we denote an ordered set of interparticle coordinate exponents, (n 0 , n 1 , n 2 , n 3 , n 4 ) which are conventionally restricted by a shell parameter Ω, If a symmetry restriction is imposed, the set of allowed {n} is constrained even more for specific values of nonlinear parameters. By construction these functions depend on twoelectron coordinates and account for the electronic correlation via explicit dependence on the coordinate r 12 . We emphasize that a construction of the basis according to Eq.
(2) allows for a very compact parametrization of the basis, because it is completely specified by just four nonlinear parameters and an integer value of Ω. To recapitulate, the trial wavefunction is represented as whereŜ ± AB = 1 ± P AB and P AB permutes the nuclei A and B,Ŝ ± 12 = 1 ± P 12 and P 12 interchange the two electrons, and appropriate ± signs are chosen to fulfil the symmetry criteria for gerade/ungerade and singlet/triplet states. By solving the secular equation one obtains linear coefficients c {n} . Such a form of wavefunction expansion is commonly referred to as the Ko los-Wolniewicz basis [12]. This nomenclature originated from the series of pioneering works of Ko los, Wolniewicz, and co-workers. It was introduced as a flexible generalization of the James-Coolidge basis aiming to represent the electron density asymmetry between the two nuclei ubiquitous in excited states of H 2 .
In previous calculations [11,12,22,23] involving tens or hundreds of KW functions in the basis, the set of {n} was carefully optimized by incremental selection of configurations that give the most significant contribution to the energy. Here we resort to a much simpler rule, given by Eq. (2), and introduce double basis functions with common values of nonlinear parameters. For convenience we introduce a shorthand notation for our basis S(Ω A , Ω B ), where S denotes a shorthand symmetry of the basis: either JC for generalized James-Coolidge (x = y = 0, u = w) [24] or KW for a general Kolos-Wolniewicz basis with no additional restrictions. Values of Ω A and Ω B define the size of sectors, which are the basis subsets carrying common independent nonlinear parameters and are constructed according to Eq. (2).
Construction of the basis by Eq. (2) for all the states considered in this work entails a universal and practical approach. It is expected, however, that more sophisticated selection of basis functions, such as subdivision into three Ω, each restricting the maximal value of exponents of η i , ξ i , and r 12 individually, as investigated by Sims [25] in James-Coolidge basis or utilized in our approach to long-range exchange splitting [26], will lead to more compact wave function expansions.

Angular factors
Expansion of the type in Eq. (3) is valid only for the states of Σ + symmetry. Extension of the KW basis to states of higher angular momentum requires the introduction of angular factors. The electronic state of a diatomic molecule is, among other symmetries, characterized by Λ -an absolute value of the eigenvalue of the n · L operator, where n is a normalized vector parallel to the internuclear axis and L is the electronic angular momentum operator. In order to enforce proper angular symmetry for Λ = 0 we construct symmetric, traceless l-th rank tensors, where ρ i = r i − n i n r, δ ij ⊥ = δ ij − n i n j , n = R/R, and (ijk) denotes symmetrization of indices. Such tensors represents irreducible representations of SO(2) rotations (around the internuclear axis) in the coordinate space. Consequently, the total wavefunction for Π, ∆, and Φ states, corresponding to Λ = 1, 2, and 3, respectively, is expanded as In the evaluation of matrix elements with the above functions, repeated Cartesian indices in bra and ket are summed over, and the resulting expression is a linear combination of f -integrals with various sets of {n}, which are defined as Efficient recursive evaluation of these integrals is the subject of our previous works [27][28][29], and the details will not be repeated here. The obtained expressions for matrix elements are too lengthy to be reported here, even those for the Hamiltonian with wavefunctions of the Π symmetry. Nevertheless, we have explicitly checked that our construction of angular factors yields the same linear combinations of f -integrals, as with the use of either real or complex angular factors introduced by Jeziorski et al. in Ref. [13]. In practical calculations, due to the length of these formulas, evaluation of matrix elements dominates the computation time for ∆ and Φ symmetry, even at moderate basis sizes.
In pioneering calculations of the lowest gerade ∆ states, namely J and S 1 ∆ g and j and s 3 ∆ g , by Ko los and Rychlewski [8,10] and later by Wolniewicz [11,14], the ππ terms were not included (the χ 12 angular factor), as hinted by Jeziorski et al. [13]. Its inclusion improves the adiabatic energies by a few cm −1 for J 1 ∆ g and S 1 ∆ g states and tenths of cm −1 for triplet j 3 ∆ g and s 3 ∆ g states, as demonstrated in subsequent work by Wolniewicz [14]. Together with the evaluation of nonadiabatic coupling by Yu and Dressler [30], the apparent discrepancy with experimental values for dissociation energies of H 2 and D 2 [31,32] has been resolved.

Results
We have calculated BO energies for 1 − 4 Π and 1 − 3 ∆ and 1 − 2 Φ states, and additionally 3 1 Φ u , and 3 3 Φ u states (38 in total) in 81 points of a nonuniform grid spanned within R ∈ (0.01, 20.0) au. In order to test the capability of the exponential basis, we have selected the I 1 Π g , i 3 Π g , C 1 Π u , and c 3 Π u states (n = 1), for which we have ventured to exceptionally large bases consisting of two sectors with their own nonlinear parameters: JC (21,19) (54648 basis functions in total) in the range R = 0.01 − 6.0 au, and KW(18,16) (53998 basis functions in total) in the range R = 6.5 − 20.0 au. This allowed us to achieve 15 or more significant digits in the range R = 0.7 − 20 au, and can be regarded as the state-of-the-art of our current computational approach.
BO energy of the 1-4 Π states In the remaining cases, calculations were performed with smaller bases, obtained by incrementing Ω of both sectors by 1 until the extrapolated energies reached uncertainty no worse than 10 −9 at all distances. This usually required the use of the JC (14,12) basis in the range 0.01 − 6.0 au and the KW(12,10) basis in 6.5 − 20.0 au. To reduce the computational cost of calculations with ∆ and Φ symmetries, both sectors having different angular factor were set to carry the same nonlinear parameters. This significantly reduced the number of f -integrals required for calculations of matrix elements.
All the potential curves have a well-pronounced single minimum around R = 2 au and exhibit strong features of anticrossings and curve interactions in the region R = 4 − 8 au, which are associated with the configuration mixing. Those interactions have a trend to become weaker with the increasing angular momentum of the excited electron. Due to the high accuracy of the calculated potential, the dominating electronic configuration can be deduced from comparison to energies of helium atom excited states and of an infinitely separated hydrogen H(1s) -H(nl) atoms. In particular, we observe, by reaching internuclear distances as small as R = 0.01 au, that the energies smoothly evolve to the corresponding value of the excited He atom [33]. It follows from analysis of curves that each highly excited molecular state is dominated by a single atomic 1snl configuration at both the united atom and the dissociation limit. The nl configuration at both of those limits is usually not the same; therefore, in such cases, when considering fixed n in the molecular term (i.e., the n-th eigenvalue of the BO molecular Hamiltonian), a character change of the excited electron has to occur at least once along the PEC. Understanding of the electronic character of a molecule in terms of singly excited 1snl atomic configurations originates from the pioneering works on molecular binding theory of Hund, Mulliken and many others [34][35][36] and can be qualitatively described with the correlation diagrams [37]. Detailed analysis of electronic characters along the PECs is out of scope of the current report and in this work we mainly focus on the computational aspects and accuracy comparison with previous ab-initio results. For a meticulous analysis of electronic configuration evolution as a function of n and R, we refer to the extensive analysis of Corongiu and Clementi [38].

States of Π symmetry
Although the basis for the Π state could in principle be built with a single sector, use of the second sector with different nonlinear parameters has proven to be advantageous. Primarily, it clearly makes the basis more flexible with the introduction of another set of nonlinear parameters. Secondarily, multi-sector bases are routinely utilized in atomic calculations with Hylleraas bases, where the largest sector spans the most diffuse length scales with nonlinear parameters bounded from below by √ −2 m E, and subsequent sectors aim to improve the shorter scale range to the nucleus. Optimal values of parameters obtained in the present calculations for H 2 suggest that spanning more diffuse functions is usually energetically more important than the short-range motion in the vicinity of the nuclei, and the introduction of a second sector seems especially relevant in the regions without the domination of a single electronic configuration.
Our numerical results for 1 − 4 Π states are plotted in Fig. 1 and presented in the Supplementary Materials in Tables S1-S16. In order to compare with the previous results of Refs. [11,15,21,39] we select R = 2 au and present value of potentials in Table 1. To test capabilities of our method, we have also evaluated a few Π u states with n = 5, 6, 10 at a single point R = 2.0 au.
In Ref. [39] an alternative approach to solving the Schrödinger equation in terms of stochastic evaluation of the path integral with the Schrödinger potential was introduced as an energy estimator. In general these results were in good agreement with previous results. However, few σ discrepancy was observed for the I 1 Π g state, for which their initial result of -0.659 515 9(3) appeared to be discrepant by a few σ with the variational result −0.659 515 055 5 of Wolniewicz [11]. Curiously enough, by taking four times larger statistical sample, more specifically four times more integration paths in the generalized Feynman-Kac method (GFK), they tightened the result to -0.659 515 54(6), thus suggesting a significant improvement over the result of Wolniewicz. Comparison of both results with ours reveals that result of Ref. [39] overshoots the variational limit by almost the same amount as the difference between our result and that of Wolniewicz. This indicates that the convergence of the GFK method to the exact energy is much slower than anticipated and the extrapolation uncertainty was much too optimistic.
In Fig. 2, we have plotted the difference between our results and that of Refs. [20,21,41,43,44] as a function of R for a few selected Π states differing in the singlet/triplet and gerade/ungerade symmetry, as well as in the character of the excited electron. Considering the energy differences between our results and those of Refs. [11,15,21] for various Π states as a function of R, we conclude that rather old results of Wolniewicz [41] for 1 − 4 1 Π u states are accurate to at least 0.1 cm −1 . Interestingly, the agreement with Ref. [21] within ∼ 1 cm −1 accuracy is observed for the 1, 2 Π states of p and d character. In contrary the discrepancy is much bigger for the states dominated by the configurations with higher l, where the discrepancy reaches even hundreds of cm −1 , see Fig 2c. Most likely it is due to the poor representation of singly-excited configurations with large values of l, especially of f, g, and h character, and could be associated with the absence of terms proportional to H(1s) -H(n = 4) in the initial function of their method. This absence is already noted by the Authors of Ref. [21].  Tables S17-S28, and are shown in Fig. 3. Moreover, in Table 2. we compare our results for ∆ symmetry at R = 2.0 au (vicinity of the minimum) with all the ab-initio results available in the literature. Comparison reveals that Wolniewicz's results for J and S 1 ∆ g and j and s 3 ∆ g states are accurate up to 0.1 cm −1 , which is roughly the uncertainty he had originally assigned to his calculations [14]. This Table reveals that for the S state, the basis JC (14,12) even with different parameters for both sectors, gives rather slow convergence to the CBS limit; therefore, the JC basis is far from optimal for this state. By virtue of high accuracy of our curves we can unambigously assign state configurations by direct comparison with helium values [33]. In particular 1, 2, 3 ∆ g and 1, 2, 3 ∆ u states start at R = 0 as 1s3d, 1s4d, 1s5d, 1s4f , 1s5f , and 1s6f , respectively and this configurations propagate at least up to R = 5 au, with the exception of 3 ∆ g states which sharply switch to 1s5g configuration at R ∼ 0.5 au. R (a.u.)

States of ∆ symmetry
BO energy of the 1-3 ∆ states   R (a.u.)

States of Φ symmetry
BO energy of the 1-3 Φ states  Tables S29-S38, and are shown in Fig. 4. In the literature Φ states were calculated only in Ref. [20] using the Free Complement local Schrödinger equation (FC-LSE) method of Nakashima and Nakatsuji [18] and in Ref. [38] using the Full-CI method. In this latter work the results were presented only for a few selected states and at specific points corresponding to the energy minimum and dissociation limit, R = 2 au and R = 100 au, respectively. Consequently, in Table 3 we present a comparison to those results at the curve minimum R = 2 au. It appears, that for the Φ state, FC-LSE functions do not span the complete Hilbert space with their choice of initial spatial functions. The deviation is slightly smaller than 1 cm −1 , and thus larger than the uncertainty estimates. Ultimately though, direct comparison with experimental data is cumbersome, because for such highly excited states the Rydberg electron decouples from the nuclear axis due to nonadiabatic effects and Λ is no longer a good quantum number. Due to high angular momentum, the curves are very regular. The difference between neighbouring states propagates almost unchanged from those of helium values [33] at R = 0 up to around R = 5 au. The singly-excited characters of 1, 2 Φ g and 1, 2, 3 ∆ u states start at R = 0 as 1s5g, 1s6g, 1s4f , 1s5f , and 1s6f , respectively and this configurations propagate at least up to R = 5 au, with the exception of 3 Φ u states which sharply switch to 1s6h configuration at R ∼ 0.3 au.

Summary and Conclusions
In order to compare our results with the plethora of transitions measured with accuracy reaching ∼ 0.001 cm −1 , theoretical values of relativistic, QED, adiabatic, and nonadiabatic corrections are necessary. The importance of the latter is arguably the most significant due to strong nonadiabatic couplings, because in contrast to the wellisolated X 1 Σ + g ground state, energy differences of BO energies between neighboring excited states can be as small as a few cm −1 . Even though some of the states with high Λ are hardly accessible with spectroscopic methods, knowledge about them has proven to be beneficial for extrapolation of high-n molecular Rydberg states. Incorporation of clamped-nuclei ab-initio potentials for relatively low n has been found fruitful in approaches based on Multichannel Quantum Defect Theory (MQDT) for Rydberg states of H 2 [45][46][47][48][49][50].
In recent years explicitly correlated Gaussian (ECG) functions have achieved a great success in high-precision calculations of hydrogen molecule isotopologues, ranging in applications from BO energy [15,51] and nonadiabatic corrections [52,53] to relativistic [54,55] and QED [56,57] corrections. All those calculations, however, were limited to the electronic X 1 Σ + g ground state. In the present work we have demonstrated that the Ko los-Wolniewicz bases still find a niche in high-accuracy variational calculations of excited states of a diatomic twoelectron molecule, especially for the states with high angular momentum. Alongside with our previous calculations of Σ + states [58], present calculations reconcile with all the previous ab-initio calculations of BO potentials of H 2 by elucidating the actual accuracy of the former calculations. Numerical uncertainty of Born-Oppenheimer potentials has been alleviated from a level of a few to the millionth parts of cm −1 ; hence, we have laid the foundation for the evaluation of further corrections to the energy levels, lack of knowledge of which currently hinders comparison with state-of-the-art spectroscopic results.

Disclosure statement
We declare no conflict of interest.