Benchmark calculations of ro-vibrational spectrum of HeH − and its isotopologues

Benchmark variational calculations of the Born–Oppenheimer potential energy curve (PEC) performed with explicitly correlated all-electron Gaussian functions with shifted centers are presented. The PEC energies include the leading relativistic and quantum-electrodynamics corrections and the first-order corrections due to adiabatic effects. The PEC is used to calculate the ro-vibrational spectra for HeH− and its isotopologues. The results show that these systems are marginally stable and have two to four bound vibrational levels and, for each vibrational level, a few bound rotational levels lying below the dissociation threshold. This indicates a possibility of detecting the HeH− anion in the laboratory and, perhaps, even in the interstellar space.


Introduction
The age of the Universe is estimated for about 13.82 billion years. In the initial stage of its existence the Universe consisted of energy in the form of photons and some elementary particles (including electrons). The atomic nuclei appeared 'only' 13.8 s after the Big Bang and they were mostly the helium, tritium, and lithium nuclei with some small amounts of the beryllium and boron nuclei. The first chemical bond that appeared in early Universe was the bond between the helium atom and a proton in the HeH + ion. This system, which is the smallest stable molecular heteronuclear ion, was recently discovered in the BGC 7027 planetary nebula located three thousand light years away from Earth in the direction of the Swan * Author to whom any correspondence should be addressed.
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. constellation [1]. As the Universe cooled down, the hydrogen atoms started to form dihydrogen molecules in the process involving interaction of hydrogen atoms with HeH + . It is the molecular hydrogen that was later responsible for the formation of the first stars. The reported observations [1] of the rotational ground-state transitions at the wavelength of 149.1 μm employed a high-altitude radio-telescope and the terahertz spectroscopy. Popular experimental methods such as isotope-ratio mass spectrometry [2] or laser-absorption spectroscopy [3,4] techniques may be useful in the HeH − anion researches.
Isotopologue abundancies are a valuable study tool for astrochemists due to the fact that they provide useful information for disentangling different chemical routes that might correspond to different physical interstellar conditions. In particular, the abundancy hydrogen/deuteron rations of hydrogen containing molecules may provide an important information of the evolution of these molecules in the Universe. It may also aid the characterization of the physical and chemical properties of these systems in the interstellar space. As the replacement of hydrogen by deuterium may induce noticeable changes in molecular geometrical parameters, chemical shifts, reaction rate constants, etc., the isotope ratio is a property worth of the investigating. Also, as the deuterium found in space was formed only during the Big Bang and in the very early Universe, and it is easy to be destroy in nuclear reactions, the ratio of deuterated and non-deuterated molecules might be a useful property that can help in exploring the origin of the Universe. The great number of deuterated molecules have been detected in recent years in the interstellar space. Theoretical calculations of shifts in spectral transitions due to isotopic substitutions may help in new detections and in the identification of new interstellar chemistries.
In view of the recent discovery of the HeH + ion [1] in the interstellar space, the existence of the negative anion formed by a helium atom and a hydrogen atom in this medium and in the laboratory can become a subject of investigation. The HeH − ion is isoelectronic with the He 2 dimer, whose existence in the phase has been experimentally confirmed. HeH − can be formed by a direct attachment of a helium atoms to a H − anion. As this heteronuclear anion has asymmetry of both its positive and negative charges, it has a non-zero dipole moment and, thus, its rovibrational spectrum should be visible in the experiment. The asymmetry can also enable the anion to radiate off the excess energy in the process of its formation from helium atoms and H − anions. Till now, the HeH − anion has been studied primarily in the context of intermolecular collisions. Most investigations performed with the use of quantum mechanical calculations were focused on the mechanism of the excess-electron removal from this system and its subsequent dissociation [5][6][7][8]. More recent studies have concerned structures of ionic clusters formed by attachment of noble gas atoms to stable positive and negative ions of the first-and secondraw elements. HeH − appears in those studies as the smallest anionic systems of that kind. See, for example, the work of Casalegno et al performed with the use of the Monte Carlo approach [9].
The bonding in HeH − is due to the electrostatic interaction between a charged atomic system, H − , and a polarizable neutral atom, He, and, thus, the main contribution to the bonding effect is the charge induce-dipole interaction. As such, the interaction is inversely proportional to the fourth power of the distance between the He and H nuclei, R −4 . As the polarizability of the helium atom is very small, the He-H − interaction is very weak and this makes the HeH − potentialenergy curve very flat and anharmonic. Naturally, this affects the ro-vibrational structure of HeH − . The depth of the potential energy curve (PEC), which is dependent on the induced dipole moment of the helium atom, can only support a very few bound vibrational levels, if any. As both, the He atom and the H − anion, are closed-shell systems, no conventional chemical bond can be formed between them and, thus, the HeH − anion can be characterized as a weakly bonded van der Waals complex. In such a complex, the equilibrium internuclear distance is usually very long and the dissociation energy is very small. Another characteristic of the complex, that can be noted, is that, while the electron density of the helium atom is compact, the H − electron density is much more diffuse. This leads to an asymmetry of the negative charge in the complex.
Before any spectral measurement is performed on the HeH − anion, it is helpful to characterize the system using quantum mechanical calculations. The first such calculations, where the HeH − PEC was determined, was the work of Olson and Liu [10,11]. Their study was done in the context of molecular collision calculations. They employed Slater orbitals and the self-consistent field (SCF) and configuration interaction (CI) methods to calculate the ground-state energy and the corresponding wave function for each PEC point. The first analysis of the stability of the isotopologues of HeH − , i.e. 4 HeH − , 4 HeD − , 3 HeH − , and 3 HeD − , was performed by Bendezzoli et al [12]. The HeH − PEC in their calculations was obtained using the full configuration interaction (FCI) method and various large Gaussian basis sets. The PEC energies were also calculated using the SCF method, the second order of Møller-Plesset perturbation theory (MP2), some modifications of the CI method, and the coupled cluster (CC) method. The calculations show that the results strongly depend on the basis set used. A more accurate HeH − PEC was obtained in more recent quadratic configuration interaction (QCISD(T)) calculations by Shalabi et al [13]. Due to inability to obtain accurately enough results for large internuclear distances, they used extrapolation to calculate the energies for these distances. The most recent calculations of the HeH − PEC are due to the work of Vallet et al [14]. Basis sets of uncontracted Gaussian orbitals were used and the results of the calculations confirmed the existence of a single bound vibrational state for 4 HeH − and 4 HeD − . Finally, in the context of the above-described HeH − calculations, one should also mention the work of Li and Lin [15], who applied a three-particle Born-Oppenheimer (BO) model to describe HeH − and to calculate the bonding energies for the different HeH − isotopologues.
Even though, the results of the above-mentioned previous calculations showed some quantitative differences, they all agreed in showing that there is at least one bound vibrational state for each HeH − isotopologue. They also agreed in describing HeH − as a system with a very shallow and anharmonic PEC and a small binding energy.

Non-relativistic BO calculations
The first step of the calculations performed in the present work involves solving the non-relativistic Schrödinger equation for every point of the PEC of the HeH − anion in the ground electronic state using the variational method and expanding the wave function in terms of, so called, explicitly correlated basis functions dependent on all inter-electron distances. The explicit dependence of the basis functions on the inter-electron distances makes it possible to achieve much higher accuracy in the present calculations than in the calculations performed before where molecular orbitals were used to construct the wave function. This is because, as the electrons repel each other, the probability of finding any two electrons close to each other is reduced in comparison with the probability of finding them more distant from each other. For example, in the two-electron helium atom, finding both electrons on one side of the nucleus is much smaller than finding the electrons on opposite sides of the nucleus. If one uses an orbital approximation and represents the electronic wave function of the atom as a product of orbitals, as it is done in the SCF method, the electron-separation effect, usually referred to as the electron correlation, is not described because the probability distribution for first electron is independent of the probability distribution of the second electron. The most effective way of making the probability distribution of the second electron in the atom be dependent on the instantaneous position of the first electron, is to include terms in the wave function that explicitly depend on the distance between the two electrons.
The BO approximation is used in the present approach and for each PEC point the corresponding electronic energy and the electronic wave function is determined. The explicitly correlated basis functions used in the present calculations are the following explicitly correlated Gaussians (ECG) with shifted centers: where A k is a positive-definite real symmetric square matrix with the dimensions 3n e × 3n e (where n e is the number of electrons) and s k is a 3n e -long vector of the shifts of the Gaussian centers. The Cartesian coordinates of the electrons form the 3n e -dimensional vector r given as: A k can be written as a Kronecker product of n e × n edimensional symmetric matrix A k with a unit 3 × 3 matrix, where ⊗ denotes the Kronecker product. A k in the above form is rotationally invariant. To ensure square integrability of function φ k (r), matrix A k must also be positive definite. This happens automatically if A k is represented in the following Cholesky factored form as: where L k is an n e × n e , rank n e , lower triangular matrix. φ k (r) is automatically square-integrable for the L k matrix elements being any real numbers. Matrix elements of L k and the elements of the shift vector, s k , are variationally optimized in the calculation for each PEC point by the minimization of the corresponding total electronic HeH − energy. In the optimization, the analytical derivatives of the energy determined with respect to these elements are calculated at each optimization step and employed by the optimization subroutine to accelerate the search for an energy minimum. The use of the derivatives, which form the energy gradient vector, is key in obtaining very accurate results, because the parameters can be more effectively optimized if the direction of an energy-lowering change of the parameters provided by the gradient vector is made know to the optimization subroutine.
The proper electronic permutational symmetry consistent with the permutational symmetry of the state of the system considered in the calculation has to be implemented. This is done through an appropriate symmetry projection operator represented by a linear combination of operators involving permutations of the electron indices. Let P be the permutation matrix with dimension 3n e × 3n e that represents permutation operator P. The following convention is used to represent the kl matrix element of an arbitrary operator O: By acting with the permutation operator, P, on basis function φ l (r) we get: whereÃ l = P T A l P ands l = P −1 s l . The procedure for deriving permutation-symmetry operators based on the Young diagrams (tableau) was described previously [16].

Relativistic corrections
Very accurate calculations of bound energy states of atoms and molecules, and the corresponding interstate transitions are possible due to very accurate values of most physical constants and due to accounting for the leading adiabatic, non-adiabatic, relativistic, and quantum-electrodynamics effects in the calculations. In the case of small atoms and small molecules formed by those atoms, these affects can be determined using the perturbation-theory approach developed by Caswell and Lepage [17]. The approach is based on expanding the energy of the system in a series of corrections proportional to positive integer powers of the fine-structure constant, α (α = 137.035 999 139): where E nr is the non-relativistic energy obtained by solving the time-independent Schrödinger equation. The α 2 correction, E rel , is an expectation value of the Breit-Pauli Hamiltonian, H rel , calculated using the non-relativistic wave function. The α 3 correction, E qed , accounts for the leading quantum thermodynamic effects. E qed , like E rel , is also calculated as a first-order perturbation-theory correction. Both E rel and E qed are corrections to the non-relativistic electronic energy, E nr . As the HeH − anion is a closed shell system, its total electronic spin is zero. Thus, the Breit-Pauli Hamiltonian that depends on the coordinates of the electrons contains the following terms: where p i is the linear momentum operator of electron i and r i is the position operator of the ith electron. H rel contains the following terms: n e j>i ∇ 2 r i 1 r i j , which represents the one-and two-particle Darwin correction, and which represents the orbit-orbit interaction correction.
The algorithms for calculating the mass-velocity correction represented by H mv and the one-and two-particle Darwin correction represented by H d were described in our previous paper [18]. Due to the fact that the previously derived algorithm for calculating the mass-velocity correction showed some oscillatory behavior for longer internuclear distances, a new algorithm is derived in the present work that replaces the old one. Due to changing the order of arithmetic operations in the new algorithm, the value of the correction now behaves smoothly with the internuclear distance. The modification, which is described below, is implemented in the computer code used in the present work. The algorithm for calculating the orbit-orbit interaction represented by H oo was presented in our previous work [19]. The algorithm is used in the present calculations.

Modification of the mass-velocity correction.
The matrix element of the mass-velocity correction is calculated as: where J ij = J ij ⊗ I 3 and J i j is a n e × n e square matrix defined as: . E i j is the n e × n e square matrix with 1 in the i j-th position and 0's elsewhere. It easy to see that: A derivation of the algorithm for calculating the mass-velocity correction in the basis set of the explicitly correlated Gaussian functions with shifted centers was shown in [18]. In order to obtain a more numerically stable value of the correction, particularly for longer internuclear distances, the new algorithm is derived for calculating the mass-velocity matrix element involving two Gaussians, φ k (r) and φ l (r). In the algorithm there is no breakdown into terms with and without the shifts (the breaking down was the source of the numerical instability of the previous algorithm). The derivation proceeds as follows: In the calculation of the above matrix element in the previous work [19] five different integrals were needed. In the present algorithm only the following two integrals are needed: and where s is a vector with length n e defined as: and Tr [A] denotes the trace of matrix A. For brevity, we use the following notation: This concludes the derivation of the algorithm.

QED correction
The total one-and two-electron QED correction of the order of α 3 for small multi-electron atoms and molecules represented by Hamiltonian H qed consists of the following terms: where • δ 3 (r i ) and δ 3 (r i j ) are one-and two-electron Dirac delta functions, Araki-Sucher term [20,21]). A matrix element of this term is calculated in the basis set of functions φ i (r) as the following integral: where Θ denotes the Heaviside step function and γ is the Euler constant. This term is calculated using the algorithm given in our previous work [22], • lnk 0 is the Bethe logarithm. For the ground state of a fewelectron atom or a molecule the Bethe logarithm can be written as: where is the ground-state wave function, N is the number of the basis functions, and c i are the linear expansion coefficients.
An algorithm for calculating the one-and two-electron QED correction of the order of α 3 (12) is implemented in the present work to preform calculations for a four-electron diatomic system with its wave function expanded in terms of all-electron explicitly correlated Gaussians with shifted centers.

Calculation of the Bethe logarithm.
The relativistic corrections (6) have been calculated for many molecular systems. However, calculations of the corrections due to the QED effects are sparse. Among the leading QED corrections, the most difficult to calculate is the term containing the Bethe logarithm. So far, this term has been only calculated for hydrogen molecular ions, HD + and H + 2 [23][24][25][26] and, in an approximate way, for the hydrogen molecule [27].
Recently, a relatively fast and effective method for calculating the Bethe logarithm, lnk 0 , was proposed by Stanke et al [28] and tested in the hydrogen-atom calculations. The method is a generalization of the method of Drake and Goldman [29], which, in addition to the method of Schwartz [30], is one of the most popular methods for calculating this term.
The starting point for developing the algorithm to calculate the lnk 0 expectation value is the theorem for calculating a function of an operator. Let us assume that the eigenvalue We assume that the solution of the Schrödinger equation H|Ψ 0 = E|Ψ 0 , is known, where HH = H 0 + H and in H and the BO approximation is assumed. Operator H is treated as a perturbation that provides a small contribution to the energy of the system. H 0 is the unperturbed Hamiltonian (here it is the non-relativistic electronic Hamiltonian). We also assume the eigenvalue problem for H 0 , H 0 |Ψ 0 = E 0 |Ψ 0 , is solved. The ground-state energy of the perturbed Hamiltonian is now written as: E = E 0 + ΔE, where ΔE = Ψ 0 |H |Ψ 0 . The expression that appears in the numerator of the expression for lnk 0 is now represented as a function of Now determination the expectation value of the above operator function with the Ψ 0 wave function is carried out. Let us take the basis set used to expand the ground-state wave function, Ψ 0 , and use it to solve the secular matrix equation for Hamiltonian H 0 : where S and H 0 are the overlap and Hamiltonian matrices, respectively, E 0 is the diagonal matrix of eigenvalues of H 0 , C is the matrix of the linear expansion coefficients of the eigenfunctions, and S = CC † −1 . Equation (16) allows to write the function of operator H 0 in a matrix form as: where f(E 0 ) is a diagonal matrix with the diagonal elements equal to the energies obtained from solving the eigenvalue equation for Hamiltonian H 0 .
To derive an expression for the Bethe logarithm we use spectral identity I constructed using an appropriate set of basis functions. The selection of the basis functions used to construct I depends on the spatial symmetry of Ψ 0 . The spectral identity, I, is inserted at two places in the expression for function f (H 0 ) of operator H 0 : From (18) it is clear that I has to be constructed using excitedstate wave functions orthogonal to Ψ 0 that have non-zero coupling matrix elements with Ψ 0 through operator ∇ r . For a diatomic molecule with the C ∞ symmetry and with Ψ 0 of the Σ symmetry, excited-state wave functions with Σ and Π symmetries have non-zero ∇ r matrix elements with Ψ 0 . Thus, in this case, the spectral identity is a sum of the Σ and Π identities: I = I Σ + I Π .
According to the procedure introduced in reference [28], the excited state basis functions can be obtained from basis functions {ϕ}, which are employed to expand the Ψ 0 wave function of the ground state, by multiplication by the appropriate spatial coordinates. In the case of a diatomic molecule with the bond axis aligned with the z-axis, two types of excite-state functions can be generated: φ Σ ∈ {zϕ} and φ Π ∈ {xϕ + yϕ}. The numbers of functions φ Σ and φ Π are n Σ = n · n e − 1 and n Π = 2n · n e , respectively, where n is the number of basis functions used to expand Ψ 0 and n e is the number of electrons in the system.
An extension of the algorithm of Stanke et al [28] to calculate the Bethe-logarithm QED term for multi-electron atoms and diatomic molecules was described in our recent work [31]. It was shown that, with the use of the extended approach, one can achieve an accuracy of two-four significant figures in calculating the Bethe-logarithm value depending on the complexity of the system. For HeH − the estimated accuracy is close to three significant figures. When the value of the correction that includes the Bethe logarithm is multiplied by the third power of the fine-structure parameter, α 3 , and added to the total energy, seven decimal figures after the decimal point in the energy value can still be considered as accurate. This allows for a reliable calculation of all ro-vibrational levels of the studied HeH − isotopologues, though, some of the very closely separated rotational levels cannot be sharply distinguished. We should note that the above-described procedure for calculating the Bethe logarithm is low-cost and can be applied in calculations for molecules larger and more complex than the hydrogen molecule. The present calculations are the first where the Bethe logarithm is determined for a molecule with four electrons using a basis set of all-electron explicitly correlated Gaussian functions.
The calculations of the Bethe logarithms for all PES points of the HeH − anion give the same value of 4.38 a.u., when the accuracy of three decimal digits is assumed. However, as the corresponding QED correction to the energy that includes the logarithm also depends of the nucleus-electron and electron-electron Dirac delta functions, the correction shows some variability with the length of the internuclear distance.

Results
As mentioned, the first calculations involve variational optimization of the explicitly correlated Gaussian basis functions with shifted centers and the calculation of the non-relativistic energy and the corresponding wave function for all PEC points of HeH − . For each PEC point the exponential parameter of the Gaussians (1), i.e. the L k and s k matrix elements, are subject to optimization. The optimization starts from the PEC point corresponding to the longest considered internuclear distance and the initial guess for the basis set for this point is constructed as a superposition of the ECG basis sets previously optimized for the He atom and for the hydrogen anion. The number of functions in the initial set is 1000. Similar approach was used by Bendezzoli et al in their HeH − calculations [12].
At large internuclear distances, the electronic wave function is less complex than for shorter distances because the helium atom can be expected to be only very slightly polarized by the distant H − . However, as the internuclear distance decreases, more electronic density starts to appear in the space between the nuclei and the wave function become more complex. This requires that additional functions are included in the basis set to maintain a comparable level of accuracy in the calculations for all PEC points. Thus, in the present calculations, the number of Gaussians in the basis set increases with the decrease of the internuclear distance. The number of ECGs at the equilibrium distance is increased to 1250. We should note that attempts are made to include additional functions for some selected PEC points corresponding to large distances (R > 50), but these additions (numbering between 10 to 100 basis functions) do not change the energies of these point by more than 10 −9 . Therefore, we conclude that the accuracy of the PEC is uniform throughout its whole length and increasing the number of the basis functions with decreasing internuclear distance does not affect the accuracy of the calculations.
After the non-relativistic energy and the corresponding wave function is determined for each PEC point, the leading α 2 relativistic (6) and α 3 QED (12) corrections are calculated. The adiabatic correction is also determined. This correction is calculated as a sum of corrections due to each of the two nuclei of HeH − according to the algorithm described in the previous paper [32]. The algorithm is based on the method developed by Cencek and Kutzelnigg [33]. Different values of the adiabatic correction are obtained for the different HeH − isotopologues.
The results for some selected PEC points obtained in the present calculations are presented in table 1. The total nonrelativistic electronic energy, E nr , and the total energy, E arq , that includes E nr , the leading relativistic and QED corrections, and the adiabatic correction for each PEC point are shown in the table. Our results obtained for the 4 HeH − isotopologue are compared with the results of Bendezzoli et al [12] and the results of Shalabi et al [13]. In general, our E nr energy values are significantly lower than the values of Bendezzoli et al and Shalabi et al. The difference appears in the third significant figure. The difference is due to much higher accuracy of the present non-relativistic energies obtained using large very well optimizes all-electron explicitly correlated Gaussian functions with shifted centers and due to including the relativistic, QED, and adiabatic corrections. In table 1 we do not show the HeH − PEC obtained by Olson and Liu [11] using the SCF method and Slater orbitals, as this PEC is much less accurate than the other PECs because it does not include the electron correlation effects. Besides the differences in the energy values between our results and the results of the others, the three PECs also differ in terms of their shapes. For example, while the equilibrium distance in our ground-state BO PEC curve of HeH − is at 12.3 a.u., the values given by Bendezzoli et al [12] obtained using the SCF, MP2, CI-SD, CC-SD, and FCI methods range between 11.41 a.u and 13.26 a.u. The SCF value of Li and Lin [15] is 11.5 a.u. The closest result to ours is the value of 12.571 a.u. from the recent work by Vallet et al [14] and obtained in calculations employing a large set of Gaussian-type orbitals.
In the next step of the calculations, the PEC's generated in the previous step are used to calculate the ro-vibrational spectra of the different isotopologues of HeH − . The calculations are done with the use of the PECs generated at the non-relativistic BO level (i.e. using the E nr energies for the PEC points), the non-relativistic level (i.e. using for each PEC point the energy, E a , being a sum of E nr and the adiabatic correction), the non-relativistic plus relativistic-correction level (i.e. the energy, E ar , for each PEC point is a sum of E nr and the relativistic correction), and at the level involving energies being sums of the E ar and the QED corrections, (E arq ). The discrete-variable-representation method [34] is used in the calculations of the bound vibrational states. The vibrational energies obtained in the calculations for the different HeH − isotopologues are shown in table 2. The dissociation energies, E d , calculated for the different vibrational states are also shown in the table. In the calculations, the mass of the helium atom is a sum of the mass of the nucleus of the respective helium isotope plus two electron masses and the mass of the hydrogen/deuterium/tritium anion is a sum of the mass of the nucleus of the respective hydrogen isotope plus two electron masses. Adding the electron masses to the masses of the nuclei in vibrational calculations allows to approximately account for the non-adiabatic effects, as these effects are due to the nuclei 'dragging' some fraction of the total electron density as they vibrate.
As the QED correction of the order of α 3 is calculated with a lower accuracy than the non-relativistic energy and the relativistic and adiabatic corrections (due to the slightly lower accuracy of the calculations of the term containing the Bethe logarithm; the number of significant figures in this term is limited to three), the calculations of the vibrational spectra with the PECs that contain all corrections (i.e. E arq ) are slightly less accurate than the calculations with the PECs without the QED correction. Thus, due to less accurate energy values of the vibrational energies in the former case, it is not always possible to determine the exact number of the bound vibrational states (i.e. states with the energies lower than the dissociation threshold). Therefore, in table 2 we only show energy values of states whose energies are fully reliable. It is very unlikely that the numbers of bound vibrational states for the different isotopologues would change when the QED correction is included in the PECs. The bound vibrational states that are obtained in the calculations, but for which the differences between their energies are below the precision of the calculations, are marked in the table with text too low precision (or tlp, for short).
In contrast to the previous works [9,[12][13][14], where only one bound vibrational state was predicted for the leading HeH − isotopologue formed by 4 He with the H − and for the other anions formed by 4 He or 3 He with H − or D − (or T − ) [12], the present calculations predict the existence of between two to four bound vibrational states. For the leading HeH − isotopologue, 4 HeH − , the present calculations predict three bound vibrational states. As expected, the number of bound  Figure 1 shows how the PEC for the 4 HeH − system shifts after consecutive corrections are added to the energy. As one can see, the addition of each corrections downshifts the PEC, but the shape of the PEC remains almost unchanged. This results in the vibrational levels shifting little upon addition of the consecutive corrections. Also, it should be noted that, as the figure shows, the depth of the PEC is much smaller than the magnitude of each correction. The positions of the bound vibrational energy levels relative to the PEC and relative to the dissociation threshold plotted for the PEC obtained from the non-relativistic BO calculations (i.e. the PEC points corresponding to E nr ) are shown in the middle section of the figure. As one can see, while the energy of ground-state vibrational level is located significantly below the dissociation energy, the energies of the two higher vibrational states are very close to each other and to the dissociation threshold (marked in orange). The central part of the figure is an enlargement of the rectangular area marked in gray in the left part of the figure In order to show the differences in energy between the three bound vibrational levels (v = 0, 1, and 2) and the dissociation energy level, the central part of the central part of the figure is enlarged and the enlarged plot is shown on the left side of the figure.
The vibrational wave functions obtained in the vibrational calculations are used to calculate the average values of the equilibrium internuclear distance for all bound vibrational states obtained in the calculations. The results are presented in table 3. Each average distance is not given as a single value but as range of values, as the error in determining the averages increases from 0.5 percent for the vibrational ground state to 2 percent for the second excited vibrational state. Also, in this case, the inclusion of the consecutive energy corrections in the PEC of the particular isotopologue has insignificant effect on the average internuclear distance. However, as expected, as the systems becomes vibrationally excited, the average distance considerably increases. This indicates that, indeed, the bond in HeH − is very weak and the anion is a very floppy structure.
In table 3 we compare our averaged internuclear distances with the results Bendezzoli et al [12] and Li and Lin [15]. As they only predicted existence of a single bound vibrational state, the comparison in only made for that state. As one can see, our values fall in between the results of Bendezzoli et al and Li and Lin.
Due to the floppy nature of the HeH − molecule, calculating bound ro-vibrational states is not as straightforward as for a diatomic molecule with a strong chemical bond. Thus, in the present work, the vibrational and rotational motions are separated, and the bound rotational levels for a particular vibrational state are calculated separately using the rigidrotor approximation. In the calculation, the moment of inertia of the molecule is determined using the vibrationally averaged distance determined for that particular state (taken form  table 3). The energy of the particular ro-vibrational state is the sum of the energy of the corresponding vibrational state obtained from the DRV calculation and the energy of the corresponding rotational state obtained from the rigid-rotor calculation.
In table 4 we show the energies of ro-vibrational states determined with respect to the ground ro-vibrational state (the state with v = 0 and J = 0) for each of the six isotopologues of HeH − . As for the pure vibrational states, the results that correspond to including all corrections in the PEC energies (due to a lower accuracy of the Bethe-logarithm term) are less accurate. Due to this lower accuracy some of the states have identical energies. This is noted in the table.
For most states, the changes in the relative energies due to including additional energy corrections are either very small or zero (within the number of significant figures used in table 4). Larger differences appear between the corresponding transition energies of the different isotopologues.
In table 5 we show a comparison of the dissociation energies corresponding to the two lowest pure rotational states, i.e. Table 6. Non-relativistic energy E nr , relativistic-E 2 α 2 , qed-E 3 α 3 and adiabatic E a (for each isotopes of helium and hydrogen) corrections in a function of inter-nuclear distance R. All values are in a.u.

R
E nr E a ( 3 He) · 10 −4 E a ( 4 He) · 10 −4 E a (H) · 10 −4 E a (D) · 10 −4 E a (T) · 10 −4 E 2 α 2 · 10 −4 E 3 α 3 · 10 −5 R E nr E a ( 3 He) · 10 −4 E a ( 4 He) · 10 −4 E a (H) · 10 −4 E a (D) · 10 −4 E a (T) · 10 −4 E 2 α 2 · 10 −4 E 3 α 3 · 10 −5 states (v = 0, J = 0) and (v = 0, J = 1), for all HeH − isotopologues obtained in the present work with the results obtained by others. All values are given in cm −1 . As the bond in HeH − is very week, the energies of all bound ro-vibrational states determined with respect to the corresponding ground ro-vibrational state are very small. None is higher than 1.3 cm −1 . Thus, it is very important that in the calculations the respective dissociation energies are determined very accurately. In the present calculations the BO PEC is calculated up to the internuclear distance of 600 a.u. At this distance the total energy is essentially constant. The comparison of the present results with the results of others [12,14,15] in table 5 shows quite pronounced differences. As mentioned before, the differences can be attributed to much higher accuracy of the present results. For example, the present value for the depth of the PEC well obtained at the non-relativistic level of 3.942 cm −1 is noticeably lower than the value of 3.732 cm −1 obtained by Vallet et al [14] and the value of 3.25 cm −1 obtained by Bendezzoli et al [12].

Summary
In this work we describe very accurate calculations of the rovibrational spectra of the HeH − anion and its isotopologues. The calculations are performed within the BO approximation and involve expanding the electronic wave function in terms of explicitly correlated all-electron Gaussian functions with shifted centers. The nonlinear parameters of the Gaussians are variationally optimized using a procedure that employs the analytical energy gradient determined with respect to these parameters. The non-relativistic energy for each point of the PEC is augmented with the adiabatic correction and the leading relativistic and QED corrections. The PEC is calculated for each of the six isotopologues of HeH − and the PECs are used to calculate bound ro-vibrational states of these systems.
The results confirm very weak bonding between the He atom and the H − ion in the system. Both, the very small dissociation energy and the large equilibrium-bond length, testify to the weakness of the bond. In this respect the present calculations confirm the conclusions of the previous works [11,12,14]. However, while the previous works predicted existence of only one bound vibrational state for the main 4 He 1 H − isotopologue, the present calculations predict existence of three bound vibrational states for this system. In general, there are between two and four bound vibrational states for each of the six isotopologues of HeH − . The accuracy of the calculations of rovibrational spectra of systems, like HeH − , is currently limited by the accuracy of determining the QED correction, in particular, the accuracy of the term that contains the Bethe logarithm. The actual accuracy of the calculation of the QED correction of the order of α 3 allows for estimating the accuracy of the PEC and of the energies of the ro-vibrational states. Nevertheless, the accuracy of the present calculations is sufficient to definitely confirm the stability of all HeH − isotopologues and to predict (with certainty) the number of bound ro-vibrational states for each of them. Also, the predicted ro-vibrational transition energies may enable spectroscopic measurement of these transitions in the laboratory and, perhaps, even in the interstellar medium, as it was recently done for HeH + .
to thank Professor PiotrŻuchowski for making his DVR computer code available to us to for his help in running the code. The data that support the findings of this study are available from the corresponding author upon reasonable request.

Appendix. Full results for HeH − molecule
See table 6.