Nonadiabatic rotational states of the hydrogen molecule
Received
22nd September 2017
, Accepted 24th November 2017
First published on 24th November 2017
We present a new computational method for the determination of energy levels in fourparticle systems like H_{2}, HD, and HeH^{+} using explicitly correlated exponential basis functions and analytic integration formulas. In solving the Schrödinger equation, no adiabatic separation of the nuclear and electronic degrees of freedom is introduced. We provide formulas for the coupling between the rotational and electronic angular momenta, which enable calculations of arbitrary rotationally excited energy levels. To illustrate the high numerical efficiency of the method, we present the results for various states of the hydrogen molecule. The relative accuracy to which we determined the nonrelativistic energy reached the level of 10^{−12}–10^{−13}, which corresponds to an uncertainty of 10^{−7}–10^{−8} cm^{−1}.
1 Introduction
The hydrogen molecule gives an opportunity to test the foundations of quantum chemistry, which are based on quantum electrodynamic theory. In principle, there are no limits to the theoretical precision of the determination of molecular levels, apart from the limited accuracy of fundamental constants, such as the electronproton mass ratio, the Rydberg constant, or the nuclear mean square charge radii. This opens up the possibility to determine these fundamental constants from molecular spectroscopic data, or alternatively to look for any discrepancies between theoretical predictions and experimental results to search for as yet unknown interactions.^{1} In fact, in recent years we have observed significant progress in the accuracy of molecular spectroscopy.^{2–10} In the particular case of the hydrogen molecule, contemporary measurements have reached the accuracy of 10^{−5} cm^{−1} (relative 10^{−9}) for selected transitions.^{11–13} On the theoretical side, various relativistic and quantum electrodynamic corrections have recently been calculated,^{14,15} but the principal problem up to now has been the insufficient accuracy of nonrelativistic energy levels. In a general multiparticle case, the complexity of the Schrödinger equation prevents its accurate solution and enforces approximations to be made. The most common one is the adiabatic approximation, which assumes the separation of the electronic and nuclear dynamics. Only a few attempts to solve directly, i.e. without the adiabatic approximation, the fourbody Schrödinger equation for H_{2} have been published. The first successful method was developed by Kołos and Wolniewicz over 50 years ago.^{16,17} They employed a nonadiabatic expansion of a trial wave function in products of electronic James–Coolidge basis functions^{18} and the vibrational functions of the form h_{n}(R) = R^{−3}e^{−x2/2}_{n}(x) with x = βR − R_{e}, where β and R_{e} are variational parameters, and _{n} denotes the nth Hermite polynomial. The expansion was composed of 54 electronic terms and six h_{n} functions yielding 147terms in total. Because of this relatively short expansion, the obtained nonrelativistic dissociation energy D_{0} = 36114.7 cm^{−1} differed by ca. 3 cm^{−1} from the exact value. Nevertheless, the pioneering work by Kołos and Wolniewicz has set the foundations of the theoretical techniques for accurate calculations and, regarding the then available computing capabilities, should be considered as a great success of theory. Fifteen years later, using the same type of wave functions with significantly larger (1070term) and carefully optimized expansion, a refined integration method, and much more powerful computers, Bishop and Cheung^{19} reduced the error to 0.2 cm^{−1}. Quite a different approach, based on the quantum Monte Carlo method, was presented by Traynor et al.^{20} in 1991 and improved later by Chen and Anderson.^{21} Their trial wave function was a product of four terms ψ_{i}. The first two terms were a combination of oneelectron functions centered on nuclei A and B, ψ_{i} = e^{−ariA} + e^{−ariB}. The third term was the Jastrow factor responsible for interparticle correlation and cusps and the last term accounted for nuclear vibration and was of the Gaussian form ψ_{4} = e^{−d(R−c)2}. The quantum Monte Carlo method allowed the finite basis set error to be eliminated but introduced instead a statistical (sampling) error, which in the latter calculations was of about ±0.2 cm^{−1}. A breakthrough result has been published by Kinghorn and Adamowicz.^{22,23} Using a 512term basis of explicitly correlated Gaussian functions, they have diminished the error in the nonrelativistic D_{0} to 1.7 × 10^{−3} cm^{−1}. Later on, successively improving the optimization technique and expanding the basis set size to 10000 terms, Bubin and Adamowicz^{24–26} arrived at an extremely accurate solution of the fourparticle Schrödinger equation to obtain D_{0} = 36118.79774(1)cm^{−1}.
All these calculations have been limited to the nonrotational state of the molecule (J = 0). In this work, we show how to incorporate the coupling between the rotational and electronic angular momentum in a straightforward manner and increase the accuracy of the nonrelativistic dissociation energy up to the level of 10^{−7}–10^{−8} cm^{−1} for the ground as well as for the rotationally and vibrationally excited energy levels of the electronic X^{1}Σ^{+}_{g} state. This project comprises one of the stages heading toward the prediction of the total energies of the hydrogen molecule with the accuracy of 10^{−6} cm^{−1}. The other contributions are relativistic O(α^{2}), leading quantum electrodynamics O(α^{3}), and higher order O(α^{4}) which are known only within the adiabatic approximation.^{15} The knowledge of nonadiabatic wave functions obtained here, is essential for the calculation of these contributions.
2 Method
2.1 From a general exponential to the nonadiabatic James–Coolidge basis function
The method described here is relevant to a molecule consisting of two electrons, labeled 1 and 2, and two nuclei, labeled A and B, with masses M_{A} and M_{B} and charges Z_{A} and Z_{B}. The nonrelativistic Hamiltonian for this system is 
 (1) 
We start the description of the method with a general exponential basis function of the following translationally invariant form

 (2) 
where
u_{1},
w_{1},
y,
x,
u, and
w are real numbers, whereas
k_{i} are nonnegative integers, and where

 (3) 
By setting u_{1} ≡ α, w_{1} = 0, y = 0, x = 0, and u = w = β we arrive at simplified basis functions

 (4) 
which still form a complete basis set. We call this function the nonadiabatic James–Coolidge (naJC) function due to its resemblance to the original James–Coolidge (JC) basis function used in fixednuclei calculations. The difference between our nonadiabatic function and the JC function is in the internuclear correlation factor as well as in the meaning of the
ζ and
η variables.
2.2 The variational nonadiabatic wave function for an arbitrary rotational angular momentum
The rotational angular momentum of nuclei couples to the electronic angular momentum, , and gives the total angular momentum of a molecule. For this reason, the wave function Ψ^{J,M} of a rotational level J (formally depending also on the projection of on the Z axis in the laboratory frame) must contain components that describe the electronic Σ, Π, and Δ,… states. In the following set of formulas we construct such a wave function and we introduce a necessary notation. The total wave function is a sum of the components with growing Λ—the eigenvalue of the · operator 
Ψ^{J,M} = Ψ^{J,M}_{Σ} + Ψ^{J,M}_{Π} + Ψ^{J,M}_{Δ}+…  (5) 
where 
Ψ^{J,M}_{Σ} = Y^{J}_{M}Φ^{J}_{Σ}, for J ≥ 0  (6) 

 (7) 

 (8) 
The particular form of functions in eqn (6)–(8) is convenient for the calculation of matrix elements as for example, the overlap matrix is block diagonal. In the above equations we use the following notation

 (9) 

 (11) 
and we assume the Einstein summation convention,
i.e. an implicit sum over all values of a repeated Cartesian index. The symbol
Y^{J}_{M} =
Y^{J}_{M}(
) denotes the spherical harmonic. The functions
Φ^{J}_{Λ} represent linear expansions in the abovedefined naJC basis functions
(4) 
 (12) 
for Λ = Σ, Π, Δ,…. For each pair
J and Λ, the function
Φ^{J}_{Λ} has its own set of nonlinear parameters, therefore we distinguish
Φ_{{k}} of
eqn (4) by indices
J and Λ. In the equation above, the symbol
_{12} means the electron permutation operator and the linear coefficients
c_{{k}} are determined variationally.
The nuclear rotation in the wave function Ψ^{J,M} is described by the spherical harmonics Y^{J}_{M}(), whereas the electronic angular contribution is represented in the form of the expansion (5) in Cartesian coordinates ρ^{i}. Each term of this expansion represents a function with a welldefined projection of the electronic angular momentum Λ. Moreover, the product of commutes with the total angular momentum operator so that it preserves the correct J and M quantum numbers. Finally, this expansion is complete. In practice, though, it can be cut due to the rapidly decreasing contribution from the subsequent terms.
A note concerning a linear dependence and a completeness of the basis set is in place here. To ensure the completeness, the function Ψ^{J,M}_{Δ} appears in two variants. The one in which both ρ and ρ′ point at the same electron, and the other, in which they point at different electrons. Certain combinations of these two variants are linearly dependent, which originates from the following identity

2_{1}_{2}(ρ^{i}_{1}ρ^{j}_{2})^{(2)} = _{2}^{2}(ρ^{i}_{1}ρ^{j}_{1})^{(2)} + _{1}^{2}(ρ^{i}_{2}ρ^{j}_{2})^{(2)}.  (13) 
This linear dependence can be avoided by, for example, using the second variant basis functions with at least one
k_{i} = 0, for
i = 2,…,5.
2.3 Symmetry of the wave function
The nonrelativistic Hamiltonian (1) is invariant under translation, rotation, and spatial inversion . The inversion reverses the sign of spatial coordinates of all particles leaving their spin unchanged, and the wave function is an eigenstate of with eigenvalues ±1. The wave function has also a definite symmetry with respect to the exchange of electrons—it is either symmetric or antisymmetric for the total electronic spin S = 0 or 1, respectively. In practice, this symmetry is enforced by acting on the spatial wave function with the operator, where _{12} exchanges the electron labels.
For a homonuclear molecule, additional symmetries arise. The gerade/ungerade inversion symmetry is the inversion of all electronic coordinates with respect to the geometric center of a molecule. Recalling how the inversion operator î acts on the electronic variables îζ_{i} = ζ_{i}, îη_{i} = −η_{i}, and îr_{12} = r_{12}, one finds that îΦ_{{k}}= (−1)^{k2+k3}Φ_{{k}} and hence

 (14) 
A wave function of a homonuclear molecule also has a symmetry due to the exchange of the nuclei. For a specified gerade/ungerade inversion symmetry and the total nuclear spin, only even or odd J levels are allowed depending on the statistics (boson or fermion) of the nuclei. For example, if we restrict our considerations to the electronic ground state () of H_{2}, we observe that the subsequent rotational levels assume alternate nuclear spins. As a result, the even J levels correspond to the nuclear singlet (parahydrogen) and the odd J levels to the triplet state (orthohydrogen).
2.4 Reduction of the angular factor
An important step in the analytic evaluation of the matrix elements with functions Ψ^{J,M}_{Λ} is the integration over the nuclear angular variables. In this section, we supply formulas for the reduction of the general matrix elements by performing the integration . Let us first note, that for an arbitrary scalar operator Q we have . In the simplest case of matrix elements with a scalar electronic (i.e. containing no differentiation over R) operator Q_{el} we have 
〈Ψ^{J,M}_{Σ}Q_{el}Ψ^{J,M}_{Σ}〉 = 〈Φ^{J}_{Σ}Q_{el}Φ^{J}_{Σ}〉  (15) 

〈Ψ^{J,M}_{Π}Q_{el}Ψ^{J,M}_{Π}〉 = 〈ρ^{i}Φ^{J}_{Π}Q_{el}ρ^{i}Φ^{J}_{Π}〉  (16) 

〈Ψ^{J,M}_{Δ}Q_{el}Ψ^{J,M}_{Δ}〉 = 〈(ρ^{i}ρ′^{j})^{(2)}Φ^{J}_{Δ}Q_{el}(ρ^{i}ρ′^{j})^{(2)}Φ^{J}_{Δ}〉  (17) 
and all the offdiagonal matrix elements vanish. The next set of formulas applies to the diagonal matrix elements with the nuclear kinetic energy operator 
 (18) 
namely 
〈Ψ^{J,M}_{Σ}∇_{X}^{2}Ψ^{J,M}_{Σ}〉 = 〈Φ^{J}_{Σ}∇_{X}^{2} − J(J + 1)R^{−2}Φ^{J}_{Σ}〉  (19) 

〈Ψ^{J,M}_{Π}∇_{X}^{2}Ψ^{J,M}_{Π}〉 = 〈ρ^{i}Φ^{J}_{Π}∇_{X}^{2} − [J(J + 1) − 2]R^{−2}ρ^{i}Φ^{J}_{Π}〉  (20) 

〈Ψ^{J,M}_{Δ}∇_{X}^{2}Ψ^{J,M}_{Δ}〉 = 〈(ρ^{i}ρ′^{j})^{(2)}Φ^{J}_{Δ}∇_{X}^{2} − [J(J + 1) − 6]R^{−2}(ρ^{i}ρ′^{j})^{(2)}Φ^{J}_{Δ}〉  (21) 
with X = A or B. Finally, the nondiagonal matrix elements read 
 (22) 

 (23) 
where + and − is for X = A and B, correspondingly, with = _{A} − _{B}. All the remaining matrix elements vanish, so that the overlap and Hamiltonian matrices have the following blockband structure 
 (24) 
2.5 Integrals with the exponential function
The previous section dealt with matrix elements without any reference to a specific shape of the spatial part of the basis function. Therefore, the above formulas can be utilized also with types of basis functions other than that presented in this article, e.g. with explicitly correlated Gaussian functions.^{26} The present section, in turn, is devoted to the exponential basis functions, in particular to the naJC functions of eqn (4). We start, however, with the most general integral for a fourbody system 
 (25) 
where 
 (26) 
The matrix elements of the Hamiltonian evaluated in the exponential basis (2) can be expressed by a combination of integrals belonging to the class defined in eqn (25). All of these integrals can be obtained through a recurrence relation starting from the so called master integral 
 (27) 
The analytical form of the integral (27) was obtained by Fromm and Hill^{27} in 1987. Their result, although terribly troublesome for a numerical evaluation, was a milestone in the evaluation of the fourbody exponential integrals. A special case of this integral was evaluated analytically by Remiddi,^{28} who expressed his result in terms of the logarithmic and the Euler dilogarithmic functions. In 1997 a significant simplification of the result obtained by Fromm and Hill was achieved by Harris,^{29} who managed to eliminate the original singularities and arrive at a much more computationally friendly formulation. Another significant step in this field was made in 2009 when the effective recurrence relations were discovered^{30} enabling evaluation of an arbitrary integral out of the whole class given by eqn (25).
The master integral g and its derivatives satisfy the following differential equations^{30}

 (28) 
where
a is one of the parameters
u_{1} ≡
t,
w_{1},
y,
x,
u, or
w, and where

σ = σ_{0} + t^{2}σ_{2} + t^{4}σ_{4},  (29) 
σ_{0} = w_{1}^{2}(u + w − x − y)(u − w + x − y)(u − w − x + y)(u + w + x + y) + 16(wx − uy)(ux − wy)(uw − xy), 
σ_{2} = w_{1}^{4} − 2w_{1}^{2}(u^{2} + w^{2} + x^{2} + y^{2}) + 16uwxy, 
The inhomogeneous term P_{a} is a combination of several logarithmic functions and is presented explicitly in Appendix A of ref. 31.
The most general integral of eqn (25) can be obtained by successive, multiple differentiation of the master integral g

 (30) 
Each differentiation increases by one the power of the associated variable in the preexponential factor of the integrand in
eqn (27). However, from the practical point of view, it is much more convenient to use recursion relations for increasing the powers
n_{i}. Let us briefly overview the steps leading to these recurrence relations. First, in
eqn (28) we set
a =
y and generate the pertinent inhomogeneous term
P_{y}. Next, we differentiate
Eqn (28)n_{2} times with respect to
y and then set
y = 0. We proceed analogously with the variables
x and
w_{1}, and obtain the relation connecting different integrals of the
Gclass. From this equation we extract
G(
n_{0},
n_{1},
n_{2},
n_{3},
n_{4},
n_{5}) and obtain the recurrence relation which, starting from
G(0,0,0,0,0,0), enables the integrals
G to be obtained for an arbitrary combination of nonnegative integers
n_{i}, expressed in terms of derivatives of
P_{y}. The multiple derivatives of
P_{y} are combinations of rational and logarithmic functions, and are numerically stable for
t − 2
u sufficiently far from zero. This condition can be easily satisfied and does not introduce limitations in practical calculations. Setting
w_{1},
y, and
x to zero simplifies significantly the analytic expressions for integrals in the naJC basis; in particular the
σ from
eqn (29) vanishes.
A small sample of explicit expressions for G and P_{y} is given below (for w = u). Note that the master integral (27) for the naJC basis is represented explicitly by G(0,0,0,0,0,0).

G(0,0,0,0,0,0) = P_{y}(0,0,1,0,0,0)/(16u^{4})  (31) 

G(0,0,0,0,0,1) =P_{y}(0,0,1,0,0,0)/(8u^{5}) + P_{y}(0,0,1,0,0,1)/(16u^{4})  (32) 

G(0,0,0,0,1,0) = G(0,0,0,0,0,1)  (33) 

G(0,1,0,0,0,0) = P_{y}(0,1,1,0,0,0)/(16u^{4})  (36) 
where

 (37) 

 (38) 

 (39) 
We note that, by symmetry, the integrals G with n_{2} + n_{3} odd vanish, as do P_{y} with n_{2} + n_{3} even. The procedure sketched above allows the whole Gclass of integrals to be evaluated analytically in a simple form.
3 Numerical approach
3.1 A perturbative solution of the eigenvalue problem
The solution of the Schrödinger equation in terms of the basis functions (4) is now reformulated into the generalized eigenvalue problem with the Hamiltonian and overlap matriceswhere is a vector of linear coefficients from eqn (12). For all states this equation can be solved directly, e.g. by the inverse iteration method. However, due to the large size of the matrices and for rotational states, it is more economical to apply the inverse iteration method only for the Σ component, and obtain the Π and Δ components from the standard perturbation theory. In other words, when J > 0, the wavefunction Ψ^{J,M}, as defined in Section 2.2, is composed of mutually orthogonal Λsegments. The orthogonality is manifested in the blockdiagonal structure of the overlap matrix, eqn (24). The small blockoffdiagonal terms of matrix enable rapidly converging perturbative expansion. In this section, we supply explicit formulas for the subsequent perturbational corrections.
Let us first consider the approximated energy level E^{(0)} = E_{Σ} obtained from the unperturbed wavefunction Ψ^{J,M} = Ψ^{J,M}_{Σ}. The Rayleigh–Schrödinger perturbation theory yields the second order (with respect to the power of the offdiagonal parts of the Hamiltonian) energy shift E^{(2)}_{Π}

 (41) 
where

 (42) 
The fourth order correction E^{(4)}_{Π} + E^{(4)}_{Δ} can be evaluated from

 (43) 

 (44) 
Each H_{Λ′Λ} Hamiltonian contains the m/μ factor of the order of 10^{−3}, which makes the perturbation series converge very rapidly. In particular, the E^{(6)}_{Δ} correction would be approximately 5–6 orders of magnitude smaller than E^{(4)}_{Δ}. Therefore, taking into account just the first four terms of the perturbative expansion is sufficient for our purposes

E ≈ E_{Σ} + E^{(2)}_{Π} + E^{(4)}_{Π} + E^{(4)}_{Δ}.  (45) 
A similar perturbation expansion holds for the wave function, namely

 (46) 
where

Φ^{(0)}_{Σ} is the unperturbed function  (47) 

 (48) 

 (49) 

 (50) 
3.2 Technical details
Each naJC basis function Φ_{{k}}, apart from the set of integers k_{i}, depends on two real positive parameters α_{k} and β_{k}. The set of the basis functions with a common pair α_{k} and β_{k} will be called the sector. Such a sector contains all basis functions with integer powers of R ranging from k^{min}_{0} to k^{max}_{0}. The total wave function can be composed of a number of such sectors. The optimal value k^{max}_{0} was determined through numerical experiments for each state separately. The k^{min}_{0} in turn was set to Λ. The ‘electronic’ integer parameters k_{1}…,k_{5} are used to organize the basis functions in ‘shells’. The given basis function Φ_{{k}} belongs to the shell number . To describe a sector of basis functions for given J and Λ values, we use a four parameter symbol (k^{max}_{0},Ω,α,β). Such a sector involves basis functions with the nonlinear parameters α and β, and with the integer powers k_{i} fulfilling k^{min}_{0} ≤ k_{0} ≤ k^{max}_{0} and 0 ≤ Ω_{k} ≤ Ω. By increasing Ω we can systematically add new basis functions to the expansion and observe the convergence of energy with an increasing total size of the basis set K.
To solve the eigenproblem (41) we employed the inverse iteration method, which consists of the ^{T} decomposition of the − E matrix followed by a solution of the linear equations set performed several times up to the assumed convergence. The matrix is block diagonal with blocks of the order 1 or 2, and is unit lower triangular. The workload of the decomposition step grows with the basis size like K^{3} whereas that of the remainder steps like K^{2}; therefore, for large matrices the decomposition step determines the timing of all of the computations. The linear algebra calculations, as well as the evaluation of matrix elements, were performed using extended precision arithmetics implemented with the help of the QD library.^{32} It enables nearly octuple precision (212 bit, 62 digits), which is sufficient to obtain the required accuracy of the energy levels considered in this work.
We consider H_{2} in its electronic ground state X^{1}Σ^{+}_{g} and, due to limited space for tables, we restrict the presentation of numerical results to the ground vibrational level v = 0. This restriction, however, is not related to limitations of the method.
3.3
J = 0 level
We consider here the ground rotational level J = 0, which requires no coupling to the electronic states with higher angular momentum to be involved, that is Ψ^{J,M} = Ψ^{J,M}_{Σ}. We used a twosector basis of Σfunctions: (30,Ω,α,β^{(1)}) and (30,Ω2,α,β^{(2)}). The parameters α, β^{(1)}, and β^{(2)} were optimized variationally with respect to the energy of the level separately for each Ω. The optimal parameters, the total size of the basis, and the resulting energy are listed in Table 1. Extrapolation of the energy to an infinite basis set size enables determination of the recommended energy value and its estimated numerical uncertainty. For this particular level, we assess the accuracy of the energy as 3 × 10^{−13} hartree. By subtracting the energy E_{0,0} from the exactly known sum of the energy of two hydrogen atoms, hartree, we calculated the dissociation energy D_{0,0} listed in the last column of the table. The numerical accuracy of D_{0,0} is estimated to be 3 × 10^{−8} cm^{−1}. This estimation, however, does not account for the uncertainty originating from determination of the fundamental physical constants. All calculations reported in this work were performed with the best currently available values of the protontoelectron mass ratio m_{p}/m_{e} = 1836.15267389(17) and of the Rydberg constant R_{∞} = 109737.31568508(65) cm^{−1} obtained from the 2014 CODATA compilation.^{33} The uncertainties of both physical constants limit the accuracy of our final results. On the other hand, the problem can be reversed and future highaccuracy relativistic calculations, in connection with highaccuracy measurements, can be applied to refine these physical constants. At present, however, the accuracy of the final value for the dissociation energy of the ground level of H_{2} is restricted by the lack of the relativistic nuclear recoil contribution and the limited accuracy of the leading QED correction.
Table 1 Convergence of the lowest eigenvalue E_{0,0} (in a.u.) and of the corresponding dissociation energy D_{0,0} (in cm^{−1}) for H_{2} with the basis set size K. A twosector wave function was employed: (30,Ω,19.19,β^{(1)}) and (30,Ω2,19.19,β^{(2)})
Ω

β
^{(1)}

β
^{(2)}

K

E
_{0,0}

D
_{0,0}

10 
0.9304 
2.664 
36642 
−1.16402503082208 
36118.797732723 
11 
0.953 
3.041 
53599 
−1.16402503087090 
36118.797743437 
12 
0.978 
3.45 
76601 
−1.16402503088047 
36118.797745538 
13 
1.011 
3.20 
106764 
−1.16402503088236 
36118.797745953 
14 
1.039 
2.80 
146072 
−1.16402503088287 
36118.797746064 
∞ 


∞ 
−1.1640250308831(3) 
36118.79774610(3) 
3.4
J > 0 levels
We consider here rotationally excited states of the ground vibrational level (v = 0, J = 1–9). In this case the admixture of states with nonzero electronic angular momentum has to be taken into account in forming the wave function (see Section 2.2 and 3.1). However, still a vast contribution of energy comes from the Σ wave function and the main effort in the calculations has to be focused in the convergence of the energy within the space formed by the Σ basis functions. For this purpose, in analogy with the J = 0 case described above, we composed the Σ wave function of two sectors of basis functions with a common α parameter: (30,Ω,α,β^{(1)}) and (20,Ω,α,β^{(2)}). By increasing the shell parameter Ω we determined the extrapolated energy value and its uncertainty. Sample data illustrating the energy convergence for a selection of states are given in Table 2. A general observation made from this table is that the convergence does not deteriorate significantly with the increasing angular momentum J, so that for the highest state considered (J = 9), the uncertainty is about the same as that for J = 1, amounting to 10^{−13} hartree. In each case, the attained accuracy in dissociation energy is better than 5 × 10^{−8} cm^{−1}.
Table 2 Convergence of the Σcomponent of selected eigenvalues E_{v,J} (in a.u.) with the increasing size of the basis set. Twosector wave functions have been employed: (30,Ω,α,β^{(1)}) and (20,Ω,α,β^{(2)}). K is the total size of the basis set and D_{v,J} – the dissociation energy in cm^{−1}
J = 1 
Ω

α

β
^{(1)}

β
^{(2)}

K

E
_{0,1}

D
_{0,1}

9 
16.92 
0.866 
2.183 
27650 
−1.1634851395410 
36000.30529283 
10 
16.93 
0.912 
2.487 
40950 
−1.1634851395785 
36000.30530106 
11 
16.99 
0.929 
3.064 
61152 
−1.1634851395847 
36000.30530243 
12 
16.93 
0.968 
3.265 
82600 
−1.1634851395863 
36000.30530277 
13 
16.93 
1.05 
3.4776 
117936 
−1.1634851395866 
36000.30530283 
∞ 



∞ 
−1.1634851395867(1) 
36000.30530285(2) 
J = 5 
Ω

α

β
^{(1)}

β
^{(2)}

K

E
_{0,5}

D
_{0,5}

9 
16.69 
0.8599 
2.282 
28756 
−1.1560957546635 
34378.52277078 
10 
15.62 
0.8998 
2.659 
42588 
−1.1560957547011 
34378.52277904 
11 
15.62 
0.9225 
3.026 
61152 
−1.1560957547073 
34378.52278040 
12 
15.62 
0.974 
3.300 
85904 
−1.1560957547088 
34378.52278073 
13 
15.62 
1.013 
3.350 
117936 
−1.1560957547091 
34378.52278079 
∞ 



∞ 
−1.1560957547092(1) 
34378.52278082(3) 
J = 9 
Ω

α

β
^{(1)}

β
^{(2)}

K

E
_{0,9}

D
_{0,9}

9 
16.45 
0.847 
2.265 
28756 
−1.1412332183826 
31116.57309932 
10 
16.29 
0.888 
2.641 
42588 
−1.1412332184218 
31116.57310791 
11 
15.39 
0.913 
3.041 
61152 
−1.1412332184285 
31116.57310938 
12 
15.39 
0.955 
2.750 
85904 
−1.1412332184298 
31116.57310966 
13 
15.39 
0.988 
3.425 
117936 
−1.1412332184301 
31116.57310974 
∞ 



∞ 
−1.1412332184302(1) 
31116.57310976(3) 
To evaluate the perturbational corrections E^{(2)}_{Π}, E^{(4)}_{Π}, and E^{(4)}_{Δ} we employed onesector wave functions of proper symmetry. This time, the sector formally depends on a single α and three β (β^{Σ}, β^{Π}, β^{Δ}) parameters. However, numerical experiments have shown that the optimal β^{Δ} is very close to the optimal β^{Σ}, and for convenience it was fixed at the value of the latter, β^{Δ} = β^{Σ}. Table 3 contains sample results of our convergence study of the three energy corrections computed according to the formulas presented in Section 3.1. An inspection of the last three columns of the table gives a view on the rate of convergence and the estimated uncertainties of particular corrections. It also informs how fast the particular corrections grow with increasing J.
Table 3 Convergence of the Σ, Π, and Δcomponents of selected rotational energy levels (in a.u.) with the increasing size of the basis set. Onesector wave functions have been employed: (30,Ω,α,β^{Σ},β^{Π},β^{Δ}). K is the total size of the basis set
Ω

α

β
^{Σ} = β^{Δ} 
β
^{Π}

K

E
_{Σ}

E
^{(2)}_{Π}·10^{8} 
E
^{(4)}_{Π}·10^{13} 
E
^{(4)}_{Δ}·10^{13} 
J = 1 
9 
19.84 
0.930 
0.875 
37429 
−1.163485134761 
−3.27271900 
−0.492866 
0.0 
10 
19.82 
0.973 
0.904 
55965 
−1.163485138021 
−3.27272362 
−0.492887 
0.0 
11 
19.83 
1.008 
0.940 
81144 
−1.163485138961 
−3.27272388 
−0.492892 
0.0 
12 
19.84 
1.049 
0.968 
114716 
−1.163485139347 
−3.27272431 
−0.492894 
0.0 
∞ 



∞ 
−1.1634851403(3) 
−3.272725(1) 
−0.492896(2) 
0.0 

J = 5 
9 
22.19 
0.923 
0.873 
58669 
−1.156095749888 
−48.9145866 
−104.5339 
−31.0603 
10 
22.19 
0.964 
0.893 
87342 
−1.156095753145 
−48.9146524 
−104.5382 
−31.0614 
11 
22.20 
0.999 
0.936 
126120 
−1.156095754085 
−48.9146553 
−104.5390 
−31.0616 
12 
22.20 
1.044 
0.964 
177644 
−1.156095754471 
−48.9146609 
−104.5395 
−31.0618 
∞ 



∞ 
−1.1560957548(4) 
−48.91467(1) 
−104.541(2) 
−31.0622(4) 

J = 9 
9 
23.90 
0.915 
0.860 
58669 
−1.141233213617 
−145.750039 
−836.830 
−259.851 
10 
23.91 
0.958 
0.889 
87342 
−1.141233216875 
−145.750228 
−836.864 
−259.857 
11 
23.91 
0.994 
0.926 
126120 
−1.141233217812 
−145.750240 
−836.871 
−259.857 
12 
23.92 
1.034 
0.953 
177644 
−1.141233218194 
−145.750258 
−836.875 
−259.858 
∞ 



∞ 
−1.1412332186(4) 
−145.75028(2) 
−836.880(5) 
−259.859(1) 
We would like to emphasize that the perturbational approach described in Section 3.1 is numerically, within the assumed goal of accuracy, totally equivalent to the variational one. The perturbational approach requires three decompositions of pertinent chunks (_{ΣΣ}, _{ΠΠ}, and _{ΔΔ}) of the Hamiltonian matrix, whereas in the variational approach the matrix must be decomposed as a whole. Because the decomposition effort is proportional to the cubic size of the matrix (∼K^{3}), the perturbational approach is, for large matrices, significantly more effective than the variational one. We have confronted the results obtained for D_{0,J} in both ways and obtained agreement better than 10^{−10} cm^{−1}. This numerical agreement shows also that consideration of only those three corrections (E^{(2)}_{Π},E^{(4)}_{Π},E^{(4)}_{Δ}) is totally sufficient for our purposes.
The final dissociation energies obtained for the lowest nine rotational levels J = 1,…,9 of the ground vibrational state are presented in Table 4. The total energy was composed of the E_{Σ} energy evaluated using a twosector wave function and the subsequent perturbational corrections E^{(2)}_{Π},E^{(4)}_{Π},E^{(4)}_{Δ} obtained from a onesector wave function. It can be seen that the final accuracy is determined mainly by the accuracy achieved for the E_{Σ} term. For higher J though, the uncertainty originating from E^{(2)}_{Π} becomes significant. The secondorder correction resulting from the coupling of the nuclear rotational angular momentum with the electronic Πstate is indispensable for accurate calculations, even for the J = 1 level. This correction increases with growing J proportionally to J(J + 1) and for J = 9 contributes to the dissociation energy as much as 0.3 cm^{−1}. The fourthorder Πstates correction and the Δstates correction are much smaller but grow even more rapidly (∼[J(J + 1)]^{2}) and become important when higherJ states or still higher accuracy is of interest.
Table 4 Nonadiabatic dissociation energy (D_{0,J}) of the lowest rotational energy levels of the ground vibrational state of H_{2}. The total as well as the Σ, Π, and Δcomponents are given (in cm^{−1}). For comparison, the results from the second order nonadiabatic perturbation theory (NAPT)^{34} are also given. The difference totalNAPT reflects the value of the higher order terms missing in the NAPT calculations
Component 
J = 1 
J = 2 
J = 3 
E
_{Σ}

36000.30530285(2) 
35764.40769523(2) 
35413.24498004(2) 
E
^{(2)}_{Π}

0.00718280 
0.02153694 
0.04304002 
E
^{(4)}_{Π}

0.00000001 
0.00000010 
0.00000038 
E
^{(4)}_{Δ}

0.00000000 
0.00000002 
0.00000010 
Total 
36000.31248566(2) 
35764.42923228(2) 
35413.28802054(2) 
NAPT 
36000.312413 
35764.429157 
35413.287941 

Component 
J = 4 
J = 5 
J = 6 
E
_{Σ}

34949.94357900(2) 
34378.52278082(3) 
33703.78059609(3) 
E
^{(2)}_{Π}

0.07165975(1) 
0.10735529(1) 
0.15007908(2) 
E
^{(4)}_{Π}

0.00000104 
0.00000229 
0.00000439 
E
^{(4)}_{Δ}

0.00000030 
0.00000068 
0.00000133 
Total 
34950.01524009(2) 
34378.63013908(3) 
33703.93068089(4) 
NAPT 
34950.015154 
34378.630045 
33703.930576 

Component 
J = 7 
J = 8 
J = 9 
E
_{Σ}

32931.16623829(3) 
32066.64655399(3) 
31116.57310976(3) 
E
^{(2)}_{Π}

0.19977878(3) 
0.25639928(3) 
0.31988489(4) 
E
^{(4)}_{Π}

0.00000759 
0.00001216 
0.00001837 
E
^{(4)}_{Δ}

0.00000233 
0.00000376 
0.00000570 
Total 
32931.36602699(4) 
32066.90296919(4) 
31116.89301872(5) 
NAPT 
32931.365910 
32066.902838 
31116.892871 
In Table 4 we compared the total nonadiabatic dissociation energy with the energy obtained from the second order nonadiabatic perturbation theory (NAPT).^{34} The difference between these two numbers comes from the higher order (m_{e}/M)^{3} terms not included in the perturbational calculations of ref. 34.
4 Summary and outlook
The numerical results presented in this work concern only the hydrogen molecule, but the method described is applicable to any fourparticle Coulomb system. The achieved accuracy surpasses that available to date for such systems. Until now, such an accuracy was available only for systems with three or fewer particles. Apart from the high accuracy, the main advantage of this method is in the formalism that enables straightforward calculations for nonzero rotational angular momentum. This feature opens up a window for accurate prediction of nonrelativistic energy for all bound levels in such systems.
The numerical results presented here constitute an introductory but indispensable part of a larger project aimed at predicting the energy levels of H_{2} with an accuracy of 10^{−6} cm^{−1}. This part must be followed by accurate (at least 1 ppm) calculations of the leading relativistic (∼α^{2}) and QED (∼α^{3}) corrections as well as corrections resulting from higherorder (∼α^{4} and α^{5}) contributions and other tiny effects like the finite size of the nucleus or gerade–ungerade mixing.^{35} We have recently evaluated the relativistic correction and the higherorder QED corrections in the Born–Oppenheimer regime.^{14,15} The new nonadiabatic wave functions will enable us to take into account also the finite nuclear mass effects in the corrections mentioned above.
Conflicts of interest
There are no conflicts to declare.
Acknowledgements
This research was supported by NCN Grants No. 2012/04/A/ST2/00105 (K. P.) and 2017/25/B/ST4/01024 (J. K.), as well as by a computing grant from the Poznan Supercomputing and Networking Center, and by PLGrid Infrastructure.
References
 W. Ubachs, J. Koelemeij, K. Eikema and E. Salumbides, J. Mol. Spectrosc., 2016, 320, 1–12 CrossRef CAS .
 J. Liu, E. J. Salumbides, U. Hollenstein, J. C. J. Koelemeij, K. S. E. Eikema, W. Ubachs and F. Merkt, J. Chem. Phys., 2009, 130, 174306 CrossRef PubMed .
 J. Liu, D. Sprecher, C. Jungen, W. Ubachs and F. Merkt, J. Chem. Phys., 2010, 132, 154301 CrossRef PubMed .
 D. Sprecher, J. Liu, C. Jungen, W. Ubachs and F. Merkt, J. Chem. Phys., 2010, 133, 111102 CrossRef PubMed .
 D. Sprecher, C. Jungen, W. Ubachs and F. Merkt, Faraday Discuss., 2011, 150, 51 RSC .
 S. Kassi and A. Campargue, J. Chem. Phys., 2012, 137, 234201 CrossRef PubMed .
 D. Sprecher, M. Beyer and F. Merkt, Chimia, 2013, 67, 257–261 CrossRef CAS PubMed .
 M. Niu, E. Salumbides, G. Dickenson, K. Eikema and W. Ubachs, J. Mol. Spectrosc., 2014, 300, 44–54 CrossRef CAS .
 M. L. Niu, E. J. Salumbides and W. Ubachs, J. Chem. Phys., 2015, 143, 081102 CrossRef CAS PubMed .
 P. Wcisło, I. Gordon, H. Tran, Y. Tan, S.M. Hu, A. Campargue, S. Kassi, D. Romanini, C. Hill, R. Kochanov and L. Rothman, J. Quant. Spectrosc. Radiat. Transfer, 2016, 177, 75–91 CrossRef .
 C.F. Cheng, Y. R. Sun, H. Pan, J. Wang, A.W. Liu, A. Campargue and S.M. Hu, Phys. Rev. A: At., Mol., Opt. Phys., 2012, 85, 024501 CrossRef .
 D. Mondelain, S. Kassi, T. Sala, D. Romanini, D. Gatti and A. Campargue, J. Mol. Spectrosc., 2016, 326, 5–8 CrossRef CAS .
 J. Biesheuvel, J.P. Karr, L. Hilico, K. S. E. Eikema, W. Ubachs and J. C. J. Koelemeij, Appl. Phys. B: Lasers Opt., 2016, 123, 23 CrossRef .
 M. Puchalski, J. Komasa and K. Pachucki, Phys. Rev. A, 2017, 95, 052506 CrossRef .
 M. Puchalski, J. Komasa, P. Czachorowski and K. Pachucki, Phys. Rev. Lett., 2016, 117, 263002 CrossRef PubMed .
 W. Kołos and L. Wolniewicz, Rev. Mod. Phys., 1963, 35, 473 CrossRef .
 W. Kołos and L. Wolniewicz, J. Chem. Phys., 1964, 41, 3674 CrossRef .
 H. M. James and A. S. Coolidge, J. Chem. Phys., 1933, 1, 825 CrossRef CAS .
 D. M. Bishop and L. M. Cheung, Phys. Rev. A: At., Mol., Opt. Phys., 1978, 18, 1846–1852 CrossRef CAS .
 C. A. Traynor, J. B. Anderson and B. M. Boghosian, J. Chem. Phys., 1991, 94, 3657–3664 CrossRef CAS .
 B. Chen and J. B. Anderson, J. Chem. Phys., 1995, 102, 2802–2805 CrossRef CAS .
 D. B. Kinghorn and L. Adamowicz, Phys. Rev. Lett., 1999, 83, 2541–2543 CrossRef CAS .
 D. B. Kinghorn and L. Adamowicz, J. Chem. Phys., 2000, 113, 4203–4205 CrossRef CAS .
 S. Bubin and L. Adamowicz, J. Chem. Phys., 2003, 118, 3079–3082 CrossRef CAS .
 S. Bubin, F. Leonarski, M. Stanke and L. Adamowicz, Chem. Phys. Lett., 2009, 477, 12–16 CrossRef CAS .
 J. Mitroy, S. Bubin, W. Horiuchi, Y. Suzuki, L. Adamowicz, W. Cencek, K. Szalewicz, J. Komasa, D. Blume and K. Varga, Rev. Mod. Phys., 2013, 85, 693–749 CrossRef .
 D. M. Fromm and R. N. Hill, Phys. Rev. A: At., Mol., Opt. Phys., 1987, 36, 1013–1044 CrossRef .
 E. Remiddi, Phys. Rev. A: At., Mol., Opt. Phys., 1991, 44, 5492–5502 CrossRef CAS .
 F. E. Harris, Phys. Rev. A: At., Mol., Opt. Phys., 1997, 55, 1820–1831 CrossRef CAS .
 K. Pachucki, Phys. Rev. A: At., Mol., Opt. Phys., 2009, 80, 032520 CrossRef .
 K. Pachucki, Phys. Rev. A: At., Mol., Opt. Phys., 2012, 86, 052514 CrossRef .

Y. Hida, X. S. Li and D. H. Bailey, QuadDouble Arithmetic: Algorithms, Implementation, and Application, Lbl46996, lawrence berkeley national laboratory technical report, 2000.
 P. J. Mohr, D. B. Newell and B. N. Taylor, Rev. Mod. Phys., 2016, 88, 035009 CrossRef .
 K. Pachucki and J. Komasa, J. Chem. Phys., 2015, 143, 034111 CrossRef PubMed .
 K. Pachucki and J. Komasa, Phys. Rev. A: At., Mol., Opt. Phys., 2011, 83, 042510 CrossRef .

This journal is © the Owner Societies 2018 