High-accuracy Rb$_{2}^+$ interaction potentials based on coupled cluster calculations

This work discusses a protocol for constructing highly accurate potential energy curves (PECs) for the lowest two states of Rb$_{2}^+$, i.e. $X\,{}^2{\Sigma}{_g^+}$ and $(1) {}^2\Sigma{_u^+}$, using an additivity scheme based on coupled-cluster theory. The approach exploits the findings of our previous work [J. Schnabel, L. Cheng and A. K\"ohn, J. Chem. Phys. 155, 124101 (2021)] to avoid the unphysical repulsive long-range barrier occurring for symmetric molecular ions when perturbative estimates of higher-order cluster operators are employed. Furthermore, care was taken to reproduce the physically correct exchange splitting of the $X {}^2{\Sigma}{_g^+}$ and $(1) {}^2{\Sigma}{_u^+}$ PECs. The accuracy of our computational approach is benchmarked for ionization energies of Rb and for spectroscopic constants as well as vibrational levels of the $a {}^3{\Sigma}{_u^+}$ triplet state of Rb\textsubscript{2}. We study high-level correlation contributions, high-level relativistic effects and inner-shell correlation contributions and find very good agreement with experimental reference values for the atomic ionization potential and the binding energy of Rb$_{2}$ in the $a\,{}^3{\Sigma}{_u^+}$ triplet state. Our final best estimate for the binding energy of the Rb$_{2}^+$ $X {}^2{\Sigma}{_g^+}$ state including zero-point vibrational contributions is $D_0 = 6179\,\mathrm{cm}^{-1}$ with an estimated error bound of $\mathcal{O}(\pm 30\,\mathrm{cm}^{-1})$. This value is smaller than the experimentally inferred lower bond of $D_0\ge 6307.5\,\mathrm{cm}^{-1}$ [Bellos et al., Phys. Rev. A 87, 012508 (2013)] and will require further investigation. For the $(1) {}^2{\Sigma}{_u^+}$ state a shallow potential with $D_0 = 78.4\,\mathrm{cm}^{-1}$ and an error bound of $\pm 9\,\mathrm{cm}^{-1}$ is computed.


I. INTRODUCTION
The recent experimental progress [1][2][3][4][5][6] towards entering the ultracold (T < 1 mK) quantum regime of hybrid ion-atom systems is naturally entangled with the need for accurate theoretical models based on first principles. While the cold regime (T > 1 mK) of ion-atom collisions is essentially classical, the ultracold domain allows for s-wave collisions and thus for reaching the pure quantum regime. The realization is still a non-trivial task for hybrid ion-atom systems due to more stringent temperature requirements, in particular for Rb 2 + [7]. Ultracold ion-atom systems provide a rich platform allowing for the discovery of novel phenomena and applications [2]. Among others, those may reach from precision measurements of ion-atom collision parameters and associated molecular potentials [1,8,9] to ultracold state-resolved quantum chemistry [10] and to the ultimate goal of realizing strongly coupled charge-neutral polaron systems [11][12][13]. For a comprehensive overview on both the theoretical and experimental state-of-the-art research on cold hybrid ion-atom systems see, e.g., Refs. [7,14].
The successful experimental realization of an ultracold ion-atom system has been demonstrated in Ref. [1] for Li 2 + . The idea of this class of experiments starts with implanting an ionic impurity into a Bose-Einstein condensate through a single precursor Rydberg atom followed by subsequent ionization to initiate collisions with the ionic core and the ground state host gas atoms. * schnabel@theochem.uni-stuttgart.de † koehn@theochem.uni-stuttgart.de ‡ lcheng24@jhu.edu For Rb 2 + , initial experiments date back to the 1960s to 1980s and have been performed in Rb vapor with densities of 2.69 · 10 19 atoms cm −3 [15]. Corresponding measurements range from associative photoionization [15,16] to rough estimates based on the analysis of charge exchange cross sections [17] to multiphoton ionization of Rb 2 and subsequent dissociation of dimer ions by one or more additional photons [18]. More recent experiments [1][2][3][4][5][6] used the same technique as described for Li 2 + . In this way, the diffuse transport dynamics of the impurity through the BEC and ion-atom-atom three-body recombination could be observed [4,5,19]. In Ref. [4] it was even possible to estimate the binding energies of some threshold bound states. These experiments may thus offer a way to probe chemical reaction channels at the quantum level. This so-called state-to-state chemistry will require resolving the quantized molecular energy levels both theoretically and experimentally.
Recent theoretical investigations of X 2 + systems (with X=Li, Na, K, Rb) can be found, e.g., in Refs. [20][21][22][23][24][25][26][27]. The equation-of-motion coupled-cluster method for electron attached states with single and double replacements (EA-EOM-CCSD) or even including triple replacements (EA-EOM-CCSDT) and including scalar-relativistic effects via the Douglas-Kroll-Hess method has been used for calculations of Li 2 + , Na 2 + , and K 2 + [25][26][27]. The authors reported satisfactory agreement with available experimental results. As perhaps the only work aiming at highest possible accuracy, Schmid et al. reported ab-initio calculations for the Li 2 + interaction potential based on an additivity scheme with coupled-cluster (CC) computations and large basis sets. These findings accompany the respective experimental work in Ref. [1], which were for the first time precise enough to predict reasonable bounds  based on the experimental approach reported in Ref. [33]. This requires the atomic Rb ionization energy EIP, the binding energy E b (see inset) of the v = 35 vibrational level of the a Σ 3 + u state of Rb2 (initial state for photoexcitation) as well as the photon energy hν required for the onset of autoionization (via a "superexcited" level of the 8 Σ for the Li+Li + scattering length. With the present work we attempt to provide a general route for obtaining high-accuracy interaction potentials for the alkaline X 2 + systems based on an additivity scheme using CC theory. The approach is inspired by the HEAT (High accuracy Extrapolated Ab-initio Thermochemistry) protocol [28][29][30], which allows for computing binding energies with uncertainties on the order of < 1 kJ/mol, which translates into < 84 cm −1 . By including our latest findings [31,32] to circumvent the wrong asymptotic behavior of some of the CC approximations required in this scheme, we report Rb 2 + PECs with so far unprecedented precision. The binding energy derived from the computed X Σ 2 + g potential can be compared to the experimentally derived [33] lower bound for D exp 0 . This experiment aimed at measuring an upper bound to the ionization energy of 85 Rb 2 formed via photoassociation, which also provides a lower bound to the dissociation energy of Rb 2 + . The bounds were measured by the onset of autoionization of the excited states of 85 Rb 2 below the 5S+7P atomic limit. As illustrated in Fig. 1, D exp 0 of Rb 2 + can be derived given the atomic ionization potential E IP , the binding energy E b [Rb 2 (v = 35)] of the initial state and the applied photon energy hν, which causes the onset of autoionization.
Besides computing spectroscopic constants, we provide analytical representations of our precise ab-initio potentials, which may serve as foundation to facilitate more sophisticated scattering calculations for more quantitative statements on the conditions needed to identify effects for Rb+Rb + collisions that go beyond the Langevin regime.
Conventional approaches [10,19,34] for studying these collision processes are merely based on model potentials.
This paper is organized as follows. Section II introduces several computational aspects of this work, including the additivity scheme proposed to potentially reach high accuracy. We gauge the expected accuracy and limitations of our approach in Sec. III in terms of benchmark calculations for ionization energies of Rb and spectroscopic constants as well as vibrational levels of the a Σ 3 + u triplet ground state of Rb 2 . Hereafter we investigate the prospects for obtaining highly accurate ground state PECs of the Rb 2 + system in Sec. IV. Therefore, we provide a protocol for constructing physically meaningful interaction potentials based on CC calculations. By thoroughly investigating basis set effects and high-level correlation contributions, we infer reasonable error bars of our computational approach and demonstrate that the latter are less important to closely reproduce experimental findings. We provide spectroscopic constants and explore the rovibrational structure of the X Σ

II. COMPUTATIONAL ASPECTS
High-accuracy quantum-chemical calculations of atomic and molecular energies often rely on additivity schemes [28,[35][36][37][38][39]. Our approach follows the HEAT protocol [28][29][30], where we assume that the total electronic energy E can be calculated via Herein, E ∞ HF and ∆E ∞ CCSD(T) are the Hartree-Fock energy and the CCSD(T) correlation energy, respectively, at the complete basis set (CBS) limit. Correlation effects beyond CCSD(T) are denoted as high-level correlation (HLC) contributions ∆E HLC and high-level relativistic (HLR) effects are labeled by ∆E HLR .
The Hartree-Fock and coupled-cluster calculations are performed using either the small-core effective core potential (scECP) ECP28MDF from Ref. [40], where only the 4s 2 4p 6 5s 1 electrons of Rb are treated explicitly and all the others are modelled via a scalar-relativistic pseudopotential (PP) or using the all-electron spin-free exact two-component theory in its one-electron variant (SFX2C-1e) [41,42] to treat scalar-relativistic effects. We use the restricted open-shell (ROHF) approach for the Hartree-Fock part and for generating the orbitals for the subsequent single-reference CC calculations. For the latter, we applied an unrestricted spin-orbital formalism in its singles and doubles variant augmented with a noniterative triples method based on a ROHF reference, the ROHF-CCSD(T) method [43][44][45][46].
The HLC contributions considered in this work include the full triples correction obtained as the difference between full CC singles doubles and triples (CCSDT) [47][48][49][50] and CCSD(T) results using smaller basis sets (i.e., e.g., triple-and quadruple-ζ quality) and the quadruples correction obtained as the difference between CC singles doubles triples augmented with noniterative quadruples [CCSDT(Q)] [51][52][53] and CCSDT results (at the triple-ζ level); cf. Sec. III. Here we used the CCSDT(Q)/B variant for ROHF references [54]. The same convergence thresholds for the HF and CC calculations as described in our previous work on the Rb 2 + system [31] were employed. The aug-cc-pCVnZ-PP and aug-cc-pwCVnZ-PP (n = 3, 4, 5) correlation-consistent basis sets for alkali and alkaline earth atoms, designed for the ECP28MDF pseudopotential [55], have been used. Furthermore, motivated by the promising results obtained in our recent work on Rb 3 [56], we constructed a large [17s14p9d7f6g5h4i] uncontracted even-tempered basis set (≡ UET17) based on the given valence basis set coming with the ECP28MDF pseudopotential in Ref. [40]; see supplementary material [57] for details. To obtain a better estimate on the Hartree-Fock limit, we combined an uncontracted and extended version of the aug-cc-pCV5Z-PP basis set for s-, p-, d-and f-exponents with diffusely augmented g-, h-, and i-functions of the UET17 basis; cf. supplementary material [57] for technical details and Tab. S.II for the respective exponents. We will call this basis "reference" in the following.
For the computations with the SFX2C-1e approach a different type of basis set is required. As the corresponding aug-cc-pwCVnZ-X2C basis sets are only available up to n = 4 [55], we used the s-, p-and d-type primitive functions of the uncontracted ANO-RCC (23s,19p,11d) set [58] and augmented them with f-, g-, h-, and i-type functions of our UET17 basis set, which were further complemented by four steeper functions to be able to describe core-valence correlation effects from the M shell (3s3p3d), see Tab. S.III of the supplementary material [57]. These basis sets will be called UANO-n with n = 3, 4, 5, 6 in this work. To finally check basis set saturation in the Hartree-Fock contribution, we replaced the s-, p-, and d-type functions by the decontracted basis functions from the aug-cc-pwCVQX-X2C basis and will call this basis "reference (ae)". The h-and i-type functions were skipped in this case to keep the size of the calculations at a manageable level. These higher-order angular momentum functions are not expected to contribute significantly at the Hartree Fock level. See Tab. S.IV of the supplementary material [57].
In detail, we used the following protocol for evaluating the theoretical best estimates of the total energy according to Eq. (1): The Hartree-Fock CBS limit is obtained by a separate computation with the 'reference' basis sets, as described above: The correlation energy at the CBS limit is computed as Here, E singles is the CCSD energy contribution resulting from the non-fulfillment of the Brillouin condition for the ROHF orbitals. We found it advantageous not to extrapolate this contribution along with the pair and triples energies as it has a distinctly different convergence behavior and saturates quickly with increasing cardinal number. Thus, we prefer to take the value for the highest available basis set result (n = 6). The pair energy (E ∞ pair ) and the noniterative triples (E ∞ (T) ) contributions are extrapolated to their respective CBS limit using the conventional two-point n −3 formula [59,60]. The same procedure was followed for ECP and SFX2c-1e computations. In the latter case, we either correlated only the 4s4p and 5s shell (the electrons explicitly treated in the ECP approach) or we also included the 3s3p3d shell (M shell) into the correlation treatment to investigate further core valence correlation effects.
The HLC contributions are computed via where we relied on ECP based calculations and the augcc-pCVnZ-PP (n = T, Q) basis sets. As indicated by the superscript "TZ/QZ", the increment for the CCSDT contribution is obtained by extrapolating the correlation energies with the aug-cc-pCVTZ-PP and aug-cc-pCVQZ-PP basis sets using the two-point n −3 formula. The CCSDT(Q) computations used the aug-cc-pCVTZ-PP basis.
The considered high-level relativistic effects are relevant as corrections for the SFX2C-1e approach only. They consist of the two-electron picture-change (2e-pc) corrections, spin-orbit (SO) corrections, and the contributions from the Breit term. In addition there are effects from quantum electrodynamics (QED), which we did not attempt to compute here. The 2e-pc correction is obtained as the difference between the spin-free Dirac-Coulomb (SFDC) [61] and SFX2C-1e results. The SO correction together with the contribution from the Breit term is calculated as the difference between the spin-orbit X2C scheme with atomic mean-field spin-orbit integrals (X2CAMF) based on the Dirac-Coulomb-Breit Hamiltonian [62,63] and the SFX2C-1e scheme. Calculations of these contributions have been carried out at the CCSD(T) level using the uncontracted ANO-RCC basis set, employing the frozen-core approximation for the 1s2s2p3s3p3d shells and deleting virtual orbitals higher than 1000 E h .
We note that it is less straightforward to incorporate ∆E HLR into ECP calculations. This is due to the fact that the ECP28MDF pseudopotential already contains a twocomponent spin-orbit coupled part with corresponding parameters adjusted to valence energies obtained at allelectron multiconfiguration Dirac-Coulomb-Hartree-Fock (DC-HF) level of theory, which includes relativistic effects at a four-component level of theory [40]. Therefore, one cannot use the difference between SFX2C-1e and highlevel relativistic effects and ECP results as the correction to ECP. Furthermore, the quality of the basis sets used in the SFX2C-1e calculations are not exactly the same as those for the ECP calculations, so taking the difference introduces additional errors.

III. BENCHMARK CALCULATIONS
In the following we benchmark the protocol described in Sec. II to assess both its expected accuracy and its limitations based on experimental reference values for Rb and Rb 2 .

A. Rb ionization potential
Calculations of the ionization potential (IP) of atomic Rb are listed in Tab. I. Comparison to the results using the aug-cc-pCVnZ-PP basis shows that the UET17 indeed leads to a much tighter convergence of the correlation energy. The gap between the n = 6 result and the estimated CBS limit of the IP is only ∼10 cm −1 and this observation also holds for the all-electron calculations with the UANO basis. We may therefore estimate the uncertainty of this result to be at most half of this difference. It is also evident from the two last columns of Tab. I that there is a significant contribution from the M shell correlation of nearly 80 cm −1 . Interestingly, the ECP result (that only implicitly includes the effect of inner shells) is rather close to the SFX2C-1e result with explicit M shell correlation, see line 'CBS' in Tab. I. As the effect of core-correlation will be much smaller in the other cases where we mainly look at binding energies, we did not go into further detail regarding this observation. Clearly, this points at the limits of the accuracy of the ECP approximation.
Moving to a larger basis for the Hartree-Fock contribution leads to corrections of the IP by 6 and 2 cm −1 for the ECP approach and the all-electron calculations, respectively. As a conservative estimate of the residual uncertainty for the Hartree-Fock contribution we thus take 5 cm −1 . The high-level correlation effects have been evaluated for the ECP/aug-cc-pCVnZ-PP computations only. The overall result is a contribution of 10 cm −1 . For the atoms, it was also feasible to carry out the CCSDT(Q) calculations with a quadruple-ζ quality basis set, which allowed a basis set extrapolation for this quantity as well, leading to another 3 cm −1 . Overall, the contribution of ∆E HLC is not very large and we will in the following use its total value as an estimate of the overall uncertainty in the correlation energy due to higher-order correlation (for the case of the IP of Rb: 10 cm −1 ). Finally, we consider the high-level relativistic (HLR) corrections beyond SFX2C-1e, ∆E HLR . The two-electron picture change correction gives 5.76 cm −1 and the inclusion of spin-orbit coupling and the Breit term contributes 2.01 cm −1 further. In addition, a QED correction of 9.98 cm −1 has been taken from Ref. [76]. These three contributions add up to ∆E HLR = 17.75 cm −1 as given in Tab. I. The uncertainty in this correction may be conservatively estimated as roughly half its amount (10 cm −1 ); we will later find much smaller relativistic corrections in the calculation of binding energies, which will then allow us to gauge the uncertainty by the total size of the high-level relativistic correction.
The final best estimates may be compared to the experimental value of E exp IP = 33690.81 ± 0.01 cm −1 [74,75]. We find a deceptively good coincidence of the ECP based result, in particular, if the extrapolated CCSDT(Q) results are taken into account. However, the error estimates discussed so far add up to approximately ±30 cm −1 . Using this error bound, the overall X2C-based theoretical best estimate of 33709 cm −1 is also in good agreement with the experimental value. The core-correlation from the M shell seems to be important for quantitative predictions and will be closely analyzed for its impact on the binding energies. It is most likely dominated by the 3d subshell, which is energetically significantly higher than the 3s and 3p subshells. The correlation from lower-lying shells (K and L) is much less important, as exploratory computations show (although with not fully sufficient basis sets, as even steeper polarization function were needed). Using the UANO-6 basis, correlation of the L shell reduces the IP by 3 cm −1 and of the K shell by further 1 cm −1 , see supplementary material [57].
In summary, the results for the ionization energy are very promising in terms of obtaining first-principle predictions with an accuracy definitely better than 1 kJ/mol (84 cm −1 ).
Next, we gauge the accuracy of our computational approach for molecular calculations and compute the spectroscopic constants D e and R e , i.e. the depth of the interaction potential and its equilibrium distance, for the lowest triplet state a Σ 3 + u of Rb 2 . This state plays a fundamental role in photoassociation processes to produce ultracold Rb 2 molecules [33,[77][78][79][80][81] and is quite challenging for computational studies due to the shallow nature of the potential, as displayed in Fig. 1. Table II lists the results in a scheme analogous to that used for the IP of atomic Rb.
Interestingly, in comparison to the known experimental reference, the ECP computations deliver the closest result. This close coincidence is largely lucky in view of the estimates for the uncertainty of these values. The uncertainty in the extrapolation is rather small in this case, the results for the CBS limit deviate by less than  0.5 cm −1 from the n = 6 value in all cases and the impact on the bond distance is less than 0.005 Å. The effect of core correlation is also very small, the results of the two sets of X2C-1e computations without and with correlation of the M shell differ by only 1 cm −1 . A slightly larger correction is imposed by the Hartree-Fock correction, which reduces the binding energy by nearly 3 cm −1 for the ECP calculations and 7 cm −1 for the all-electron calculations. The high-level corrections for full three-electron clusters and four-electron clusters work in opposite directions, giving corrections of −2 cm −1 and +6 cm −1 , respectively. High-level relativistic corrections turn out to have very little impact on the binding energy, they reduce the binding energy by 0.30 cm −1 and thus do not close the small gap between the ECP and the all-electron results. In view of the smallness of this correction (a similar size of correction will be found for Rb 2 + ) we decided to not include this value in the total sum for the best estimate but only as indicator for the size of effects due to the approximate treatment of relativistic effects. The total estimate of the uncertainty of the binding energy, obtained in the same way as before for the IP, adds up to ±12 cm −1 and thus also spans the difference between the ECP and all-electron calculations. The experimental value lies quite at the border of this range, indicating that the all-electron calculations may have the general tendency of overbinding due to basis set superposition errors. The predicted equilibrium bond distance follows the trend for the binding energies: The larger the binding energy the shorter the bond. Again in comparison to the experimental value, the all-electron results give slightly too short bond distances.
We also assessed the size of the spin-orbit splitting of the triplet state by performing frozen-core MRCI/ECP28MDF/UET17(n = 3) calculations [83][84][85] using the ECP-LS technique for the corresponding smallcore ECP. The computations included seven singlet and seven triplet states of Rb 2 , such that in total a 28×28 spin-orbit matrix is set up and diagonalized. The resulting zero-field splittings and energy shifts for the a Σ 3 + u state are both 1 cm −1 , which suggests that spin-orbit effects can be neglected in the further discussion. A very small zero-field splitting could also be observed in the X2CAMF computations. This is in accordance with the results reported in Ref. [86] for the a Σ 3 + u state of Rb 2 , yielding effects in the order of 0.5 cm −1 for D e when comparing numbers obtained with and without SO effects included.
The ROHF-CCSD(T)/ECP28MDF/UET17 approach was also used to compute a larger set of points of the potential energy curve and the one-dimensional RP-RKHS interpolation method as described in Refs. [87][88][89][90][91][92][93] was used to obtain an analytic representation. The respective longrange coefficients required by this procedure were taken from Ref. [80], with C 6 = 0.227003 · 10 8 cm −1 Å 6 , C 8 = 0.778289 · 10 9 cm −1 Å 8 and C 10 = 0.286887 · 10 11 cm −1 Å 10 ; cf. Sec. IV B for details on the asymptotic form of interaction potentials. This long-range behavior with three reciprocal power terms defines the RP-RKHS parameters, yielding n Rb2 = 3, m Rb2 = 2 and s Rb2 = 2 as well as R a = 18.0 a 0 . To further provide a physically meaningful short-range description, we imposed a correction of the form This only affects the PEC for R < 4.0 Å and is introduced to correct for inherent artefacts of the RP-RKHS procedure concerning short-range extrapolation [94]. Thus, it has no impact on D e and R e but gets relevant for high-energy scattering states (E > 15000 cm −1 ). Inspired by the approach reported in Ref. [93], we further provide an empirically adjusted "optimal" PEC for which we scaled and shifted the CBS ab-initio data in order to match the experimentally derived values D exp e and R exp e , respectively [32]. The resulting PEC can be reproduced using our programs and data available from Ref. [95]. In the present work, we used this PEC to compute the rovibrational structure of the a Σ 3 + u state of Rb 2 . The resulting 41 vibrational levels for J = 0 are given in Tab. S.VIII of the supplementary material [57]. We observe excellent agreement with experimentally measured levels.

IV. HIGH-ACCURACY RUBIDIUM ION-ATOM INTERACTION POTENTIALS
The previous discussion on Rb ionization energies and spectroscopic constants of the a Σ 3 + u triplet ground state of Rb 2 demonstrate the capability of our computational approach to predict energies and potential energy curves to high accuracy. In the following we investigate the binding energies and the full PECs of the two lowest states of Rb 2 + .
A. Accurate binding energies of Rb2 + We start out with a discussion of the binding energy and equilibrium distance of the X Σ 2 + g state. As already discussed previously [31], there is in our approach a slight consistency error in the asymptotic region, as the symmetry-adapted mean-field (Hartree-Fock) solution formally dissociates the system into two fragments which are a 50:50 mixture of a cation and a neutral atom. Thus, the computed energies in the asymptote do not coincide with the sum of those of the atom and the cation from individual mean-field computations. The consistency could only be achieved by allowing for symmetry-broken solutions in the asymptote. This problem may be fully avoided by special techniques, e.g. the electron attachment equation of motion coupled-cluster (EA-EOM-CC) approach [25][26][27]. However, going to higher-order correlation methods, the effect diminishes and eventually does not significantly contribute to the overall uncertainty of the final result. At the Hartree-Fock level, using the ECP based approach, the asymptote lies 40 cm −1 above the atomic limit, which diminishes to −4.4 cm −1 at the CCSD(T) level (nearly independent of the basis set size). At the CCSDT and CCSDT(Q) level, the error shrinks to +1.5 cm −1 and +0.7 cm −1 respectively. In the SFX2C-1e scheme, the same observations are made, differing by only 0.3 cm −1 from the values quoted above. We will thus in the following report energies relative to those computed for the separated ion and atom. The main results are summarized in Tab. III, using the same scheme as discussed before.
The CCSD(T) binding energies at the CBS limit differ by approximately 20 cm −1 when evaluated from either the ECP or the SFX2C-1e all-electron approach. When correcting for the Hartree-Fock limit, using again large basis sets in both approaches, the difference shrinks to less than 10 cm −1 . Correlating the M shell in the latter approach adds 14 cm −1 to the binding energy. The difference between ECP and SFX2C-1e results cannot be attributed to remaining high-level relativistic effects. By comparison of SFX2C-1e and SFDC computations, we can estimate a two-electron picture-change effect on the order of −0.14 cm −1 whereas going to the two-component X2CAMF scheme gives a correction of the dissociation energy by −0.26 cm −1 , resulting in a total ∆E HLR of −0.4 cm −1 . The high-level correlation effects are also rather small. Including full triply connected clusters decreases the binding energy by 3 cm −1 and the correction from CCSDT(Q) increases it again by approximately half of this value.
Overall we may estimate an uncertainty of ±10 cm −1 for the Hartree-Fock contribution, ±5 cm −1 for the extrapolation of the CCSD(T) CBS limit and ±5 cm −1 for high-level correlation effects. Further uncertainties may be correlation effects from even deeper core shells (should not exceed ±5 cm −1 ) and further relativistic effects (which are also unlikely to exceed ±5 cm −1 ). Our best estimate for the dissociation energy of Rb 2 + in the X Σ 2 + g state is thus 6202±30 cm −1 . The predicted equilibrium distance changes by less than 0.02 Å among all calculations, except for those with the smallest basis set, see Tab. III. Correlation of the M shell leads to a small contraction of the distance by 0.01 Å and our best estimate is thus 4.796(10) Å.
The antisymmetric (1) Σ 2 + u state also shows a shallow minimum, which is purely due to induced dipole-charge and dispersion effects. The first conclusion is supported by the fact that the computed Hartree-Fock contribution to the binding energy is around 60 cm −1 at the equilibrium distance. The best estimate for the binding energy including correlation effects can be deduced from the numbers in Tab. IV. In this case there are only little effects from basis set extrapolation, both for the correlation and the Hartree-Fock contribution. Likewise there is only a small discrepancy of 1 cm −1 between ECP-based and SFX2C-1e-based computations and a nearly vanishing effect of core correlation. High-level relativistic effects have not been estimated for this case due to technical problems in converging the orbitals, but based on previous experience they are expected to be miniscule. The largest correction comes from high-level correlation effects, a reduction of the binding energy by nearly 5.5 cm −1 from the CCSDT calculations and a slight increase again by 2.9 cm −1 from CCSDT(Q). The final best estimate is thus 80±9 cm −1 for the binding energy.
The equilibrium distance also shows only very little dependence on the basis set size and the biggest correction comes from the high-level correlation effects, being a shortening by 0.015 Å. There is also some effect of M core correlation, which shortens the bond length by 0.02 Å in comparison to the computation that keeps this shell uncorrelated. Overall, we may give a best theoretical estimate for R e of 12.163(20) Å.

B. General long-range form of X2 + interaction potentials
An essential aspect of the recent experimental works [1][2][3][4][5][6] towards entering the ultracold domain of ion-atom interactions relies on studying scattering events. The scattering properties of such collisions are defined by the long-range form of the respective interaction potentials [7]. For ionic dimers with a single active electron (i.e., e.g., alkali-metal systems such as Rb 2 + ) this long-range behavior contains the two contributions [14] V where V ind/disp (R) describes the leading induction and dispersion interaction, while V exch defines the exchange interaction term. For interactions between an S-state atom In leading order the interaction is due to the charge q of the ion inducing an electric dipole moment of the atom, which reflects in the corresponding induction coefficient with the static electric dipole polarizability α d of the atom. The next higher-order induction term is caused by the interaction between the charge of the ion and the induced electric quadrupole moment of the atom, described by the respective coefficient with the static electric quadrupole polarizability α q of the atom. The van der Waals type dispersion interactions are usually weaker, with the leading-order term of Eq. (7) accounting for dynamic interactions due to instantaneous dipole-induced dipole moments of the ion and the atom arising due to quantum fluctuations. Higher-order terms could be added to Eq. (7), but since reaching the s-wave scattering regime for Rb 2 + is barely possible at all, truncating after the R −6 is usually sufficient.
The second term of Eq. (6) is due to the indistinguishability of the two limiting cases Rb+Rb + and vice versa. It is thus defined by the energy splitting between the asymptotically degenerate gerade and ungerade states and determines the cross section for resonant charge transfer [96]. For alkali-metal X 2 + systems, it generally involves one active electron and if both ion and atom are in an S-state, it takes the form, in atomic units [14,20] In Ref. [14], the parameters α, β and B are related by simple expressions to the ionization potential I Rb of the Rb atom The corresponding value for the ionization potential is taken from Ref. [74,75], yielding I exp Rb = 0.15350655 E h (original measurement: I exp Rb = 33690.81 ± 0.01 cm −1 ). The parameter A is the normalization factor of the asymptotic wavefunction involved to arrive at Eq. (10b) and is given [20,97] in terms of the parameters of Eq. (11) where Γ(·) denotes the Gamma function. The secondorder expansion coefficient C of Eq. (10b) may be extracted from fits to ab-initio results and is taken from Ref. [14] with C = −19.22.
The characteristic ion-atom interaction length scale R * may be derived from Eq. (7), yielding [7] R * = 2µC 4 −2 , which is in general at least one order of magnitude larger than corresponding neutral atom-atom interactions. The corresponding characteristic energy scale E * ∝ (2µ 2 C 4 ) −1 is at least two orders of magnitude smaller than the one for neutral atom-atom systems. Due to the small reduced mass of Li 2 + , this explains why it was possible to reach the s-wave scattering regime in the experimental run described in Ref. [1]. For the Rb 2 + system, the characteristic interaction length scale is R * ≈ 5000 a 0 , with the respective s-wave scattering limit E * = k B × 79 nK [2]. This stringent temperature requirement is one of the reasons why reaching the quantum collision regime for Rb 2 + is considerably more difficult as compared to neutral atom systems or to Li 2 + .

C. Construction procedure
The results from Sec. IV A indicate that overall the ROHF-CCSD(T) approach based on effective core potentials is already quite accurate. We will thus use it for the overall construction of the X Σ 2 + g and (1) Σ 2 + u potential energy curves and account for higher-level corrections by an appropriate rescaling.
Concerning the CCSD(T) approach for X 2 + systems in general, our previous work [31] revealed some limitations, as already indicated at the start of this section. In Ref. [31] we demonstrated that CCSD(T) leads to an unphysical long-range barrier of the respective system, which is related to a symmetry instability of the underlying Hartree-Fock mean-field solution. However, our findings also suggested that using (T) corrections from symmetry-broken calculations for the long-range tail and properly merging these with symmetry-adapted solutions for smaller internuclear distances may be a promising approach to construct well-defined and physically meaningful global PECs for the X Σ with the hybrid (T) correction to model the long-range region given by where Θ(R) is the Heaviside function. The first two terms in Eq. (13) denote the ROHF reference energy and the CCSD correlation energy, respectively. The point group labels represent the computational point groups and correspond to symmetry-adapted (D 2h ) and symmetry-broken (C 2v ) solutions. Equation (14) formally represents the use of (T) corrections from symmetry-broken calculations to model the long-range tail and the proper merging to symmetry-adapted (T) corrections at some merging point R m in the intermediate region to describe the remaining part of the PEC. At R m the symmetry-adapted values have to be shifted by the respective constant energy difference |∆E s | to the symmetry-broken solution to obtain continuous curves. This construction procedure is schematically illustrated in Fig. 2. In total energies this shift amounts up to ≈ O(14 cm −1 ), as can be approximated by the respective differences in Fig. 2. For interaction energies, in particular in the vicinity of the potentials equilibrium R = 4.8 Å, the respective difference is in the order of ≈ O(0.5 cm −1 ). In general, the point R m , cf. Eq. (14), enters as an additional free parameter that might be adjusted so to optimally reproduce certain experimental findings (e.g. spectroscopic constants, scattering length, vibrational levels, etc.).  We further note, that the small oscillation that occurs in Fig. 2 for the symmetry-broken curve before it collapses to the symmetry-adapted one is a consequence of the numerical bistability of the symmetry-broken CCSD(T) solution in the region where it cannot "decide" whether to collapse to the symmetry-adapted Σ g or Σ u state; see Refs. [31,32] for details.
The X Σ 2 + g and (1) Σ 2 + u states of Rb 2 + are asymptotically degenerate and thus the construction procedure according to Eqs. (13) and (14) can be applied to obtain physically correct PECs for both states. This implies that in both cases R m has to be chosen equally, which is only justified if Σ g and Σ u are reasonably well degenerate for a given merging parameter R m . Here, we assume that this is fulfilled if the difference between the total symmetryadapted ROHF-CCSD(T) energies is ≤ 10 −8 E h . Moreover, the choice of R m should be sufficiently far off the repulsive long-range barrier occurring at R ≈ 100 Å [31]. These two conditions define the restriction We set this value to R m = 40.0 Å for the following discussion. As shown in the previous section it is important to extrapolate the ROHF-CCSD(T) results to the respective basis set limits. Again, we found an approach based on Eqs. (2) and (3) as the best compromise to obtain the CBS values using the UET17 basis sets. Effects beyond the CCSD(T)/ECP approach, as discussed in Sec. IV A, are now simply included by rescaling the potential energy curves to match D e and R e of the theoretical best estimates.
To generate PECs that can be used to study Rb + +Rb scattering events it is inevitable to recover the correct longrange behavior according to Eq. (6). This also involves ensuring that the rescaled ab-initio data reproduce the theoretically suggested exchange splitting between the gerade and ungerade states as given by Eq. (10b). The following protocol is used to incorporate the exchange interaction into the PEC construction procedure [32]: 1. Compute the difference ∆(R) between the theoretically suggested exchange splitting V exch (R) of Eq. (10b) [with corresponding parameters from Eqs. (11) and (12)] and the one resulting from the rescaled ab-initio data V exch according to Eq. (10a)

Equation (10b) is only valid for [97]
Rβ 1 and In the following we assume that this means that both products have to be at least one order of magnitude larger than one. This yields a validity interval V of internuclear distances R within which the rescaled PECs have to be modified. Note that this is independent from the merging parameter R m as introduced above to join symmetrybroken and symmetry-adapted (T) solutions in the long-range tail.
3. The long-range parts V LR of the X Σ and where R v labels the lower bound of the validity interval (18).
An overview of exchange splittings V exch that result from ab-initio data and a corresponding comparison with the respective theoretical curve according to Eq. (10b) is given in the supplementary material [57]. Moreover, it is shown that the exchange interaction is very sensitive to the respective basis set leading to an interchange of the Σ g and Σ u states; a problem that already occurs at the SCF level. Finally, the adapted ab-initio data (≡ "hybrid ROHF-CCSD(T)/CBS/mod" level of theory), i.e. those obtained by first transforming according to the hybrid model of Eqs. (13) and (14), followed by a rescaling to the theoretical best estimates for D e and R e and then by subsequent modification to obtain the correct exchange splitting using Eqs. (19) and (20), are passed to the one-dimensional RP-RKHS interpolation method [87][88][89][90][91][92][93]. This gives analytical representations for the X Σ 2 + g and (1) Σ 2 + u PECs and yields by construction the correct leading-order multipole terms according to Eq. (7). It is important to include a sufficient amount of training data in the region R ∈ [20.0, 60.0] Å to ensure precise interpolation of the exchange interaction [32]. Beyond that region, the training data should be chosen sparsely to ensure extrapolation after Eq. (7). The corresponding induction coefficients follow from the experimentally measured static electric dipole polarizability α d [98,99] yielding [see Eq. u state. Furthermore, Eq. (5) was imposed to provide a physically meaningful short-range description. Again, the RP-RKHS parameters and the shortrange correction do not affect D e and R e . The resulting PECs can be generated using the python script and data available from Ref. [101], which are either based on the aug-cc-pCVnZ-PP or the UET17 basis set families both extrapolated to their respective CBS limit. Corresponding data were generated with the hybrid ROHF-CCSD(T) approach according to Eqs. (13) and (14) and modified to reproduce the correct exchange interaction using Eqs. (19) and (20). This procedure yields the PECs depicted in Fig. 3.

D. Resulting potential energy curves
As displayed in Fig. 3, the ground state of gerade symmetry is much deeper compared to the asymptotically degenerate state of ungerade symmetry, which only shows a rather shallow well at about 12 Å. We note that the choice of the merging parameter R m , within the constraints defined by Eq. (15), mainly affects the long-range behavior of the respective PECs. It can thus be viewed as defining a lower bound of the fitting range if the corresponding ab-initio data were used to extract higher-order induction and dispersion coefficients (i.e. C disp 6 , C ind 8 , etc.). These higher-order coefficients might improve the quality of the RP-RKHS interpolated PECs by including them into the inherent extrapolation according to Eq. (7). This in turn may be used to screen the sensitivity of subsequent scattering calculations based on these RP-RKHS PECs that account for such higher-order effects. A more detailed discussion on the extraction of induction and dispersion coefficients from hybrid ROHF-CCSD(T) ab-initio data may be found in Ref. [32]. As already discussed in Sec. IV C, the choice of the merging point R m and thus of |∆E s |, cf. Eq. (14), leaves R e effectively unchanged and alters D e by about 0.5 cm −1 . [32]. The corresponding effects that result from accounting for the theoretically suggested exchange interaction are also very small and give rise to the tiny differences between the theoretical best estimate and the model potential for the binding energies D e as documented in Tab. V. The rovibrational ground state energies, i.e. (v, J) = (0, 0) where v and J denote the vibrational and rotational quantum numbers, were extracted using the Level16 code [102] and the final RP-RKHS potential energy curve. The results are summarized in Tab. V.
Experimental data given in this table refer to the works already mentioned in the introduction. The most recent experiment [33] aimed at measuring the ionization potential of 85 Rb 2 formed via photoassociation of ultracold 85 Rb atoms, cf. Sec. I. The molecules were subsequently excited by single-photon UV transitions to states above the ionization threshold. This approach yielded an upper limit for the ionization energy of 85 Rb 2 and simultaneously provided a lower bound for D 0 of 85 Rb + 2 ; as reported in Tab. V with D exp 0 ≥ 6307.5 cm −1 . As indicated in Tab. V, theoretical works on Rb 2 + found in the literature can be divided into two classes: calculations that are based on model potentials [23,103,104] or approaches using large-core effective core potentials (lcECPs) with a core polarization potential (CPP) [22,33,[105][106][107] to account for polarization effects of the core electrons. In the latter approach spin-orbit effects can be incorporated into the pseudopotential, cf., e.g., Refs. [22,106]. The full valence configuration interaction ab-initio calculations that accompany the experimental work in Ref. [33] were based on a relativistic non-empirical large-core pseudopotential with CPP. Spin-orbit effects were modeled through a semi-empirical spin-orbit pseudopotential and taken into account in the valence CI calculations via the CIPSO procedure (configuration interaction with perturbation including spin-orbit coupling) [86]. As can be inferred from the numbers reported in Ref. [22], the SO effects on the spectroscopic constants of the two lowest states of Rb 2 + are negligible.
In comparison to the most recent experimental reference value D exp 0 ≥ 6307.5 cm −1 (Ref. [33]) the binding energy computed in this work comes out too low by approximately 120 cm −1 . While this is in general a satisfactory agreement, the discrepancy is somewhat in conflict with the accuracy estimates given by us previously, which amount to O(30 cm −1 ).

E. Rovibrational structure
The analysis of the rovibrational term values supported by the X Σ g state amount to about 46 cm −1 , close to values reported from earlier computations [33].
The (1) Σ 2 + u potential is extremely shallow as reflected by ω e ≈ 3.8 cm −1 . Despite its well-depth of only ≈ 80 cm −1 it still can support more vibrational levels than the a Σ 3 u state of Rb 2 (cf. Tab. S.V) with a potential depth of ≈ 241 cm −1 , for which only 41 bound states are found. This is a direct consequence of the exceedingly large interaction length scale R * of Rb 2 + , which is at least one order of magnitude larger as compared to the neutral species [7].
All data generated by the Level16 code can be found in the supplementary material [57], which also includes the corresponding input files for technical details and to rerun our calculations. An analogous study based on the aug-cc-pCVnZ-PP basis set series and referring to Rb 2 + containing merely the 87 Rb isotope may be found in Ref. [32].

V. CONCLUSIONS
This work provides a protocol for computing highly accurate binding energies and accurate global PECs for the ground state of Rb 2 + within an additivity scheme based on ROHF-CCSD(T) calculations. The approach circumvents our recently revealed limitations of perturbative coupled-cluster approximations by using a hybrid model with symmetry-broken (T) corrections describing the respective long-range part of the PEC properly merged to symmetry-adapted solutions for smaller internuclear distances. Furthermore, the construction procedure is designed such that the PECs for the X Σ states of Rb 2 + reproduce the physically correct exchange splitting. This is particularly important when using the corresponding potentials for highly accurate studies in the context of ultracold chemistry, e.g. for scattering calculations. In this regard, we moreover provide readyto-use analytical PEC representations obtained within the framework of RP-RKHS interpolation.
We benchmarked the accuracy of our computational method for ionization energies of Rb as well as for spectroscopic constants and vibrational levels of the a Σ 3 + u triplet ground state of Rb 2 . In both cases we obtained very good agreement with respective experimental findings, which confirm the tight error bars estimated for the uncertainty of the computed values. While in this respect including high-level correlation (HLC) contributions and high-level relativistic corrections are essential for ionization energies, we found that they are of much less importance for interaction energies. For both Rb 2 and Rb 2 + , HLC effects are in the order of ≈ 2 − 3 cm −1 for the total binding energy. Our computational approach gives a value of 6202±30 cm −1 for the dissociation energy of the X Σ state. We also investigated the rovibrational structure of these states based on the model potentials fitted to our computed values using the RP-RKHS approach. In particular we get a lowest rovibrational state in the X Σ 2 + g with a binding energy of 6179 cm −1 . There is a residual deviation of approximately 130 cm −1 from the most recent experimental estimate for the lower bound on the binding energy of Rb 2 + . This is a closer agreement as compared to other theoretical works on X 2 + systems, with X ∈ {Li, Na, K, Rb}. However, the discrepancy is significantly larger than the estimate for the residual uncertainty of our computations. This will certainly require further theoretical and experimental investigations.
The data that support the findings of this study are available within the article and its supplementary material [57].