Pre-Born-Oppenheimer energies, leading-order relativistic and QED corrections for electronically excited states of molecular hydrogen

For rovibronic states corresponding to the $B$ and $B'\ ^1\Sigma_\text{u}^+$ electronic states of the hydrogen molecule, the pre-Born--Oppenheimer (four-particle) non-relativistic energy is converged to a 1-3 parts-per-billion relative precision. The four-particle non-relativistic energy is appended with leading-order relativistic, leading- and estimated higher-order quantum-electrodynamics corrections. The resulting term values referenced to the rovibronic ground state are obtained in an excellent agreement with the experimental results. Further results are reported and discussed for other rovibronic states assignable to the $C\ ^1\Pi_\text{u}$ and the $EF,GK,$ and $HH\ ^1\Sigma_\text{g}^+$ electronic states.

H(1)+H (2) H(1)+H(3) Figure 1. Rovibronic states computed in this work can be assigned to the electronic states B, B ′ 1 Σ + u , C 1 Πu, and EF, GK, HH 1 Σ + g , for which the potential energy curves are visualized in this figure [19][20][21]. The shaded background of the gerade manifold indicates that the corresponding four-particle states are (predissociative) resonances embedded in the H(1)+H(1) continuum. The T BO energy is defined with respect to the minimum of the X 1 Σ + g electronic ground state.
In this work, we report computational results for the rovibronic states corresponding to the singlet ungerade manifold with N = 0, 1, . . . , 5 total angular momentum quantum numbers, p = (−1) N natural parity, v = 0 and 1 vibrational and B electronic state labels ( Figure 1). Furthermore, computational results for the vibrational 'ground state' of the B ′ electronic state is also reported, which was obtained as an excited state within the (N, p, S e , S p ) = (0, +1, 0, 1) (total angular momentum, parity, total electronic and protonic spin quantum numbers) 'symmetry' block [12], in which the v = 0, B state is the absolute ground state. The computed four-particle non-relativistic energies are appended with perturbative relativistic and quantum electrodynamics (QED) corrections computed in this work and using the Bethe logarithm values taken from Ref. [13].
The present work is an extension and continuation of four-particle computations carried out for the v = 0, N = 0, 1, 2, . . . , 5, EF 1 Σ + g inner-well rotational states [14], obtained in an excellent agreement with experiment [6,8]. Precise computation of electronically excited rovibronic states of the singlet gerade manifold is challenging, due to coupling to the lower-energy H(1) + H(1) dissociation channel. The lowest-energy EF inner-well states could be computed almost as bound states due to their very long predissociative lifetime [15]. Initial results and further thoughts for other electronically excited singlet gerade states are reported at the end of this work, in relation with EF , GK, and (inner-well) HH states. Interestingly, the outer-well of HH 1 Σ + g is relatively weakly coupled to other nearby electronic states as it was pointed out by Wolniewicz, Ubachs, and co-workers in 1997 [1,2]. In 2019, we used single-state non-adiabatic perturbation theory to improve [16] upon the adiabatic result [1,2]. Further improvement would be possible by multi-state nonadiabatic perturbation theory [17,18] or perhaps by a highly efficient four-particle predissociative computation (e.g., by further development of Refs. [12,14]).

A. Variational solution of the four-body Schrödinger equation
We consider H 2 = {e − , e − , p + , p + } as a four-particle system, and solve the translationally invariant Schrödinger equation, with n + 1 = 4 where Hartree atomic units are used, m i and q i are the masses and electric charges of the particles, The energy (and wave function) for a selected (kth) eigenstate can be converged to high precision using a non-linear variational procedure [22][23][24]. Each eigenstate can be characterized by the total spatial angular momentum quantum number and corresponding projection, N, M N , the parity, p, and the spin quantum numbers of the electrons, S e , M Se , and of the protons, S p , M Sp . Hence, the wave function is approximated by a linear combination of anti-symmetrized (A) products of φ where λ = (N, M N , p) and ς = (S e , M Se , S p , M Sp ) collect the spatial and spin quantum numbers. The c i linear combination coefficients are computed by solving the linear eigenvalue problem over the selected basis set. The χ [ς] spin function is the product of the singlet electron and singlet or triplet proton spin functions depending on the spatial symmetry of the computed state [12,23] ( Table I). For spatial functions, we use explicitly correlated Gaussian functions with the global vector representation (ECG-GVR) [23,25] (for the precise and fully detailed definition, please consult Ref. 23), where the so-called global vector is and the spatial direction of v,v = v/|v| determines the polar angles of the Y N MN (v) spherical harmonics function. The non-linear parameters α ij , u i , and the integer value K = 0, 1, . . . are optimized and fine-tuned to minimize the (kth) non-relativistic energy. For fine-tuning the real parameters (α ij and u i ), we used the Powell method [26]. Analytic matrix elements for ECG-GVR functions and the non-relativistic operators have been implemented in a Fortran computer program in Ref. [23], and this implementation has been used in practice (in double precision arithmetic in Fortran) up to ca. K = 20 and up to N = 5 − 6 [12,14,23,24]. For the computations and the non-linear parameter optimization, the translationally invariant kinetic energy operator is expressed in terms of atoms-in-molecules internal coordinates [23].
At this point, we mention two important alternative ECG basis sets relevant for diatomic computations. First, ECG functions with r 2K 1 -type prefactors have been efficiently used to compute diatomic molecules with N = 0 [27][28][29] and N = 1 [30][31][32]. For N = 0, this basis set can also be considered as a special case of ECG-GVR with restriction of the global vector to contributions from the two protons only. Currently, our ECG-GVR implementation includes contributions from all particles, which is necessary to have a fundamentally correct total angular momentum representation beyond N = 0. Second, complex ECG (CECG) functions have been proposed [33,34] to efficiently describe the nodal structure of nuclear localization and excitation. In spite of the appealing simplicity of the CECG integrals, tight optimization of the non-relativistic energy with respect to the CECG (non-linear) parameters is an intricate task, due to numerical instabilities, yet CECG implementations exist for N = 0 (and 1) [33][34][35][36][37].

B. Leading-order relativistic and QED corrections
In the non-relativistic quantum electrodynamics approach [38,39], the energy is written as the sum of the non-relativistic energy and correction terms for higher orders of the α fine-structure constant (in this work, the α powers are counted according to the use of Hartree atomic units), For computing relativistic and QED corrections to the pre-Born-Oppenheimer (preBO) energy, we consider electronic quantities (along a series of fixed nuclear configurations). The reference state is obtained by solving the electronic Schrödinger equation (for fixed nuclear configurations, with r 12 = |r 12 | = |r 1 − r 2 |, r ia = |r ia | = |r i − R a |, and R = |R 1 − R 2 |. The φ electronic wave function is approximated as a linear combination of anti-symmetrized products of two-electron singlet spin and spatial functions. For spatial functions, we use floating explicitly correlated Gaussian functions (fECGs), which we adapt to the 1 Σ + u or 1 Σ + g symmetry [13] of the relevant electronic state considered in this work, and optimize the s i ∈ R 2×3 and A i ∈ R 2×2 (symmetric, positive definite matrix) parameters based on the variational principle.
The dominant correction due to special relativity is obtained at α 2 -order as expectation value of the Breit-Pauli Hamiltonian, which, for a singlet φ electronic state, reads as with including the so-called mass-velocity, one-electron Darwin, two-electron Darwin plus spin-spin Fermi contact interaction, and orbit-orbit correction terms, respectively. The spin-independent α 3 -order QED correction [40,41] is commonly written in a compact form as [38] ε where is the (non-relativistic) Bethe logarithm, where φ and ε is the non-relativistic wave function and energy corresponding to the H el electronic Hamiltonian, Eq.
The so-called Araki-Sucher term [40,41] (related to the retardation correction in the intermediate photon momentum range [41]) can be written in the usual form as with the Θ step function and the γ E Euler constant. Furthermore, higher-order quantum electrodynamics corrections in this non-relativistic perturbative framework can be estimated as [42,43] ε (4) hqed,est = π Then, we compute our current best energy value for the kth rovibronic state as the sum of the preBO energy, E k (which we henceforth label as E (0) k to emphasize that it is a non-relativistic quantity) from Eq. (1), and the expectation value of the ε (2..5) κ relativistic and QED corrections (for the corresponding, dominant κth electronic state) with the Ψ k preBO wave function, where r 34 highlights the proton-proton distance dependence of the correction, which is integrated with the four-particle wave function to obtain the correction. At this point, it is necessary to mention that the four-particle Breit-Pauli Hamiltonian expectation value [44] and even a four-particle Bethe logarithm value [43] has been reported for the rovibronic ground state (X 1 Σ + g ) of H 2 . The present approach is much simpler and it can be expected to give useful results, if the electronic state provides a meaningful zeroth-order approximation or in other words, if the non-adiabatic-relativistic (and QED) coupling is small. The overall accuracy of the present results is thus limited by the magnitude of this coupling.
The preBO energies were obtained and converged within repeated variational non-linearoptimization cycles of the ECG-GVR basis set, separately carried out for every rovibronic state. Since the low-energy (ro)vibronic states corresponding to B and B ′ 1 Σ + u are bound states, we could use the standard variational optimization procedure applicable for a selected (kth, discrete) state [22] (Sec. IV A). Tight variational convergence of electronically excited rovibronic states from the the singlet gerade manifold is more challenging due to predissociation, initial results and thoughts are collected in Sec. IV B.
After the variational non-linear optimization procedure of the preBO energy had been sufficiently converged (including a total of 9000 basis functions optimized separately for every state), this energy was appended with relativistic and QED corrections computed as expectation value of the ε (n) (R) (n = 2, 3, 4, 5) correction curves with the four-particle wave function, Eq. (17).
As to the leading-order relativistic electronic correction, Wolniewicz reported values for the B and B ′′ B 1 Σ + u , as well as for the EF, GK, HH 1 Σ + g electronic states in Refs. [3,52] and [53], respectively, but we could not find relativistic correction values for the B ′ 1 Σ + u state. Regarding the leading-order quantum electrodynamics correction, the i,a δ(r ia ) and δ(r 12 ) values can be calculated from the one-and two-electron Darwin terms (of the relativistic correction), but we could not find any literature data for the small, though non-negligible Araki-Sucher term for these states.
To have the corrections up to α 3 -order, we computed the Born-Oppenheimer electronic (energies and) wave functions, and all relativistic and QED terms, Eqs. (10), (12), (14)-(16) for the B and B ′ 1 Σ + u , as well as (with possible future relevance), for the EF, GK, and HH 1 Σ + g electronically excited states. To pinpoint the expectation value of the singular operators (mass-velocity, Dirac delta functionals, and Araki-Sucher term) with the electronic wave function described by the current fECG basis set (including 1200 functions), we used the integral transform (IT, regularization) technique originally proposed by Pachucki, Cencek, and Komasa [45,47,54]. The relativistic and QED correction values (obtained with the IT technique) are sufficiently accurate to obtain the term correction to 1−5·10 −9 E h precision [45,47]. All computed data is deposited in the Supplementary Material.
The computationally demanding Bethe logarithm has been recently reported by two of us for the relevant states near their equilibrium structure(s) [13], and we used the value for the dominant electronic state at its (relevant) equilibrium structure (dependence on the proton-proton distance and non-adiabatic coupling was neglected).
During this work, we used the CODATA18 recommended values for constants and conversion factors [55].

IV. NUMERICAL RESULTS
A. Selected rovibronic states corresponding to the singlet ungerade electronic manifold a. B and B ′ states The non-relativistic preBO energies and term values computed for rovibronic states assigned (based on their energy) to the B and B ′ 1 Σ + u electronic states are collected in Table I and are estimated to be converged within 5 · 10 −9 E h (0.001 1 cm −1 ) for the v = 0 and within 1 · 10 −8 E h (0.002 2 cm −1 ) for the v = 1 vibrational levels. The reported experimental uncertainty for the observed term values [6] is 0.000 1-0.000 5 cm −1 . Table II collects the relativistic and quantum electrodynamics corrections to the computed term values, up to α 3 order and estimates for the α 4 -and α 5 -order contributions. The final computed energy has an estimated numerical uncertainty determined by the convergence of the preBO energy, furthermore, the overall accuracy is limited by the magnitude of the non-adiabatic-relativistic coupling (which is difficult to estimate without computation).
The overall agreement of the computed results with the experimental values, within the uncertainty of both, is excellent and represents significant improvement with respect to earlier work of Wolniewicz, Orlikowski, and Staszewska including non-adiabatic coupling as well as relativistic and estimated QED corrections [56], and the earlier, purely non-adiabatic work of Wolniewicz and Dressler [57]. We also note that Dressler and Wolniewicz [58] reported a very good estimated theoretical term energy (not shown in the table) for the B rovibrational ground state, obtained from a clever combination of various relevant corrections available at that time. It is interesting to note in Table II that the α 3 -order corrections are (in an absolute value) larger than the α 2 -order contributions to the term energy for every reported B state. This relation is reversed for the B ′ state, though the α 3 -order correction remains significant.
b. C state Finally, between the B and B ′ 1 Σ + u potential energy curves, we can observe the potential energy curve of the C 1 Π u electronic state in Figure 1. We report for the rovibronic 'groundstate' corresponding to C 1 Π u , the non-relativistic preBO energy (a variational upper bound to the exact energy) obtained and optimized as an excited state within the (N, p, S e , S p ) = (1, −1, 0, 0) 'symmetry block' to be E  Table I. Non-relativistic preBO energy, E (0) in E h , and term value, T (0) in cm −1 , for B and B ′ rovibronic states with N = 0, 1, . . . , 5 total angular momentum quantum number, natural parity, p = (−1) N , and v = 0 and 1 vibrational state labels. The electron spin is Se = 0, the proton spin is Sp = (1 + p)/2, i.e., Sp = 1 proton triplet for even N and Sp = 0 proton singlet for odd N values.  Table II for α n ε (n) X X00 (n = 2, 3)), we find that their non-relativistic energy value (not reported separately in Ref. [56]) is by 0.028 cm −1 larger, than our variational result. Computation of the relativistic and radiative corrections to our variational four-particle energy is left for future work, and it will allow direct comparison with the experimental result [6], which is currently by ∆T (0) C00 = −0.143 cm −1 lower, than our non-relativistic term energy. Table II. Perturbative relativistic and QED corrections computed to the non-relativistic preBO term energies (Table I) We report initial results and observations regarding selected rovibronic states from the singlet gerade (EF, GK, HH) manifold. The electronic expectation values for the relativistic and QED corrections (most importantly, the Araki-Sucher term) are deposited in the Supplementary Material, and the Bethe logarithm for these states near the equilibrium structures are available from Ref. [13]. For computation of the preBO energy, an initial basis parameter set was compiled from earlier bound-state computations [12,14,23,24]. Then, the basis representation of the (long-lived, but predissociative) rovibronic states was optimized in a stabilization-like procedure [14].
a. Inner-well HH state Regarding the inner well of the HH 1 Σ + g electronic state, we could converge the preBO energy for the rotational-vibrational 'ground' state (v = 0) as excited state within the (N, p, S e , S p ) = (0, +1, 0, 0) symmetry block, it is E  Table I). This non-relativistic term value differs from experiment [6] by ∆T (0) o-c = 0.125 cm −1 , which is less than the 0.89 cm −1 deviation of the non-relativistic, non-adiabatic result reported by Yu and Dressler [60]. For the relativistic and QED corrections, we could not pinpoint the numerically stable result using the procedure of Sec. II B and Eq. (17), and direct computation of the four-particle relativistic corrections appears to be necessary.
b. GK states Regarding the optimization of the lowest-energy vibrational states corresponding to the GK 1 Σ + g electronic state, we could observe the 'G0' (v = 1 for GK) and also the 'K0' (v = 0 for GK) states, but their non-relativistic energy could be obtained only with a ca. 1 µE h uncertainty in a stabilization-like preBO computation [14]. Thus, our current best estimated four-particle energies are E  [6,61], which is smaller than the deviation of the non-adiabatic value, 3.38 and 2.635 cm −1 , reported by Yu and Dressler [60]. Further improvement of the optimization procedure, perhaps in combination with optimization of the complex-energy of the complex-coordinate rotated Hamiltonian [23,62] could be used to pinpoint the four-particle energy (as a resonance state).
c. EF states As to rovibronic states corresponding to the EF 1 Σ + g electronic state, rotational excitations of the vibrational ground state (E0N ) of the 'E' inner well have been reported to high precision including relativistic, QED, and higher-order QED corrections in an excellent agreement [14] with the (more accurate) experimental results. In this work, we could also converge the preBO energy of the first vibrationally excited state assigned to the inner well (E10, but v = 3 for EF ), it is E o-c,E10 = 0.188 cm −1 experiment-theory deviation, which is smaller than the 0.746 cm −1 value from the non-adiabatic computations of Yu and Dressler [60]. In further work, relativistic and QED corrections to this state should be obtained in direct, four-particle computations.
Regarding the 'F ' outer-well states, our current best four-particle energy for the rovibrational ground state (F 00, but v = 1 for EF ) is E F 00 = 99 365(1) cm −1 , which differs from the experimental value by −0.7 cm −1 . This deviation is larger than the 0.14 cm −1 difference of the non-adiabatic result of Yu and Dressler [60]. In our procedure, the non-linear optimization of this state is stable, but currently inefficient. To speed up the computation, it will be necessary to find a better starting basis set (the current starting parameterization was compiled from basis functions optimized for rovibronic states with ca. 2 bohr proton-proton expectation value, which is much smaller than the proton-proton distance characteristic for F states).

V. SUMMARY AND CONCLUSION
Rovibronic states from the singlet ungerade (B, B ′ , C) and the singlet gerade (EF, GK, HH) manifolds of molecular hydrogen have been computed by solving the electron-proton four-particle Schrödinger equation. For the lowest-energy rovibronic B and B ′ states, an excellent agreement of the rovibronic term value (referenced to the absolute ground state) is obtained with respect to high-resolution experiment after inclusion of the relativistic, leading-order and estimated higherorder quantum-electrodynamics corrections. For the other states, the current status and necessary further developments are discussed.