Quantum Data-Driven Modeling of Interactions and Vibrational Spectral Bands in Cationic Light Noble-Gas Hydrides: [He2H]+ and [Ne2H]+

1. Introduction

Molecular spectroscopy is a fundamental tool for investigating molecular structures, properties, dynamics, and chemical reactivity. However, spectral data do not directly reveal molecular information; instead, it is encoded in complex patterns that require careful analysis. To extract meaningful insights, computational modeling plays an essential role in decoding these spectra. By combining advanced methods and techniques, such as quantum chemistry calculations, molecular dynamics simulations, and machine learning algorithms, computational spectroscopy allows researchers to simulate and interpret spectral features, enabling a deeper understanding of molecular behavior. Recent advancements in computational methods have significantly expanded the applicability of such approaches, facilitating in this way direct comparisons with experimental data [1,2,3,4,5,6,7,8]. A computational model is built upon specific concepts and assumptions, ideally adaptable to any molecular system. Quantum mechanical approaches, in particular, rely on the molecular Hamiltonian operator, incorporating both kinetic and potential energy, to extract vibrational information. Achieving high accuracy in spectral simulations requires both rigorous nuclear quantum simulations and high-level electronic structure calculations to provide an accurate description of interaction energies and an efficient representation of the potential energy surface (PES). The PES is a cornerstone for understanding most molecular processes, and its precise representation is crucial for reliable spectral computations. Commonly, discrepancies between theoretical predictions and experimental observations arise from inaccuracies in the underlying potential representation. As molecular complexity increases, refining PES models becomes progressively more challenging.

This study focuses on computing molecular interactions and characterizing the vibrational transitions of cationic light noble-gas hydrides of astrophysical interest. For a long time, the interstellar medium (ISM) has been considerate inhospitable environment to chemical compounds due to extreme conditions that characterized it and, only a few molecules could be form, but recently over 300 molecules have now been identified [9]. Notably, the discovery of noble gas hydride cations like the ${\mathrm{HeH}}^{+}$ in the planetary nebula NGC 7027 [10] and ${\mathrm{ArH}}^{+}$ for first time in Crab Nebule [11], and later, in extragalatic sources [12], challenges the notion of noble gases as inert, indicating the need of further research for a much better understanding of the noble gas chemistry.

Among the cationic noble-gas hydrides, the most astrophysically relevant are those of helium and neon. The interest in [${\mathrm{He}}_2{\mathrm{H]}}^{+}$ lies in its composition of helium and hydrogen, the two most abundant elements in the universe, increasing the likelihood of its formation in suitable space environments. The fifth most abundant element is the neon, and thus it is also an intriguing candidate. However, no molecular species containing it has been detected so far, the NASA’s LADEE mission recently confirmed [13] the presence of neon in the lunar exosphere, suggesting that neon hydrides or their cations could exist in cold planetary systems under favorable conditions for Ne interactions with H or ${\mathrm{H}}^{+}$.

In the laboratory, high-resolution mass spectrometry experiments [14,15,16,17] have provided insights into the stability of specific structures, revealing differences with those observed in heavier noble gas hydrides, [${\mathrm{Ng}}_n{\mathrm{H}}_m{\mathrm{]}}^{+}$ (m=1,2 and 3). Thus, noble gas matrix experiments has confirmed the existence of centrosymmetric [${\mathrm{NgHNg]}}^{+}$ configurations [18,19,20,21], highlighting that symmetric insertion of a proton into a weakly bound noble gas dimer significantly enhances the system’s stability compared to other possible arrangements. Just recently, applying different action spectroscopic techniques [8,22] the ${\nu }_3$ proton shuttle motion of the ${\mathrm{He}}_2{\mathrm{H}}^{+}$ has been observed at $1315.84$, together with seven of its rovibrational lines [8], while estimates for the bending and symmetric stretch fundamentals have also been reported [8,22,23] at 874.9 and $956.6$, respectively.

Several theoretical studies have been also carried out predicting both geometric and energetic characteristics of the linear configurations of these systems, identifying the centrosymmetric structure as the most stable [24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,40]. In particular, for the [${\mathrm{He}}_2{\mathrm{H]}}^{+}$ case, Dykstra [27] has first computed potential energy values for a limited set of geometries, which later on, Lee and Secrest [37] have used them to parameterize an analytical PES and calculate rotational-vibrational states up to J=2. They have identified five bound states for J=0 for both [${\mathrm{He}}_2{\mathrm{H]}}^{+}$ and [${\mathrm{He}}_2{\mathrm{D]}}^{+}$. Subsequently, Baccarelli et al. [28] have explored the potential for both collinear and ${\mathrm{C}}_2$ symmetry geometries using MRCI/pVTZ calculations, while another significant study was carried out by Kim and Lee [29], examining a restricted region of the configuration space with the CCSD(T) method. In 2003, Panda and Sathyamurthy [34] constructed a PES using the many-body expansion form [41] fitted to data from CCSD(T)/AVTZ calculations, and reported seven and fourteen bound vibrational states for the [${\mathrm{He}}_2{\mathrm{H]}}^{+}$ and [${\mathrm{He}}_2{\mathrm{D]}}^{+}$, respectively. Later, Liang et al. [38] have developed another analytical PES using 15682 MRDCI/AV5Z data points up to energies of $10000$. More recently, Stephan and Fortenberry [32] have reported quartic force fields derived from CCSD(T)/AV5Z calculations, while Fortenberry and Wiesenfeld [33] have reported an update PES using CCSD(T)-F12 data. In contrast, for the [${\mathrm{Ne}}_2{\mathrm{H]}}^{+}$ complex the literature is more limited. There are only two PESs available: a quartic force-field PES [31] using CCSD(T)/AV5Z calculations, and an analytical 3D PES [35] constructed from more than 23000 CCSD(T)/AVQZ data points. In this latter study the first 22 vibrational bound levels have been also reported for the [${\mathrm{Ne}}_2{\mathrm{H]}}^{+}$ [35].

Advances in computational quantum chemistry coupled with cutting-edge experimental techniques have opened new avenues to explore the fundamental chemical properties of these species. A common aspect in all previously mentioned computational studies is the use of traditional analytical PES models. While, such models have demonstrated the ability to provide high-quality results [42,43,44,45,46,47,48,49], even today, despite the feasibility of high-level quantum chemistry calculations, their construction still remains a laborious and time-consuming computational task. In this context, machine learning (ML) PESs trained on previously calculated molecular structures and energies [50,51,52], presents a promising alternative [53,54,55,56,57,58] for studying various physical-chemical processes, e.g. the vibrational analysis of molecules. In our recent investigations [59,60,61,62], we have assessed the performance of machine learning PESs methods based on both kernel and neural networks (NN) representations [56,63,64,65,66,67]. Specifically, for small size systems, the reproducing kernel Hilbert space (RKHS) method [68,69,70] has proven effective, as by construction accurately reproduce the training data, while simultaneously capturing potential’s long-range asymptotic interactions using appropriate kernel polynomials [71,72,73,74,75,76].

Therefore, in this work, we aim to develop a full 3-dimensional kernel-based RKHS ML-PES for the [${\mathrm{He}}_2{\mathrm{H]}}^{+}$ and [${\mathrm{Ne}}_2{\mathrm{H]}}^{+}$ noble gas hydride cations. These PES models will be trained and validated against high-accuracy CCSD(T)/CBS data, with each step of the PES construction process being thoroughly examined. The resulting ML-PES will be utilized in ro-vibrational quantum computations to provide precise predictions of energetic and spectroscopic properties for common isotopologues, which will then be compared with available literature data. The calculated vibrational transitions could provide valuable insights for supporting the astrochemical identification of these noble gas compounds in unexplored regions of the interstellar medium.

2. Computational Methods, Results and Discussion

2.1. Electronic Structure Calculations: Reference Data on Interaction Energies

All ab initio electronic structure calculations were performed using the MOLPRO 2022 program [77], while the DENEB software package [78] was employed to generate and organize all input and output data files, respectively.

2.1.1. Optimized Structures and Dissociative Energetics

We have first carried out geometry optimizations of these complexes at the CCSD(T)/aug-cc-pV6Z level of theory together with the corresponding harmonic vibrational normal modes frequency analysis. Figure 1 summarizes both energies and configurations of stationary points and dissociation channels. We have found that the optimized equilibrium structure (global minimum) of both triatomic complexes show a linear ${\mathrm{Ng-H}}^{+}$-Ng configuration, belonging to $D_{\infty h}$ symmetry, with the Ng–H⁺ bond lengths of 0.925/1.143 Å and a total energy of -5.9030/-257.8380 a.u. for the He/Ne-systems. We can note that the Ng–H⁺ bondlengths are larger compared with 0.77/0.90 Å in the isolated [${\mathrm{NgH]}}^{+}$ molecules [79] due to the incorporation of the new Ng atom. We also found that the antilinear ${\mathrm{Ng-Ng-H}}^{+}$ optimized structure (local minimum) is at energy of -5.883120/-257.8137 a.u., 4367.16/$5326.65$ above the global one, while the barrier between them is found at energy of -5.883095/-257.81372 a.u, or 4368.64/$5331.039$.

The harmonic vibrational frequencies for the ⁴${\mathrm{He-H}}^{+}{\mathrm{-}}^4$He and ²⁰${\mathrm{Ne-H}}^{+}{\mathrm{-}}^{20}$Ne from the CCSD(T)/AV6Z computations are shown in Table S1 (see Supplementary Material), where $v_1$ is the symmetric ${\mathrm{Ng-H}}^{+}$ stretch, $v_2$ the doubly degenerate proton bend vibration, and $v_3$ the asymmetric ${\mathrm{Ng-H}}^{+}$ stretching or shared-proton mode. The calculated harmonic zero-point-energy (ZPE) values is $2295.065$ for [${\mathrm{He}}_2{\mathrm{H]}}^{+}$ and $1782.385$ for the [${\mathrm{Ne}}_2{\mathrm{H]}}^{+}$, with the harmonic vibrational frequencies being in line with the reported values in previous studies [31,32,36]. Regarding the corresponding intensities, that of shared-proton stretching modes have higher values compared to $v_2$ modes, indicating a strong absorption or emission line in the IR-spectra. Furthermore, one can be seen that the proton-shuttle motion is very similar in both species, with a difference of just $3$. We also note that, in the case of [${\mathrm{He}}_2{\mathrm{H]}}^{+}$, the trend in frequencies observed in its heavier counterparts is not followed, where $v_1$ < $v_2$ < $v_3$ is the rule, instead, the order is partially inverted ($v_2$ < $v_1$ < $v_3$) and, the value of $v_1$ is twice that in the Ne-complex, and almost six times larger than in the ${\mathrm{Ar}}_2{\mathrm{H]}}^{+}$ [60].

Next, we also performed CCSD(T)/AV6Z computations to determine the energetics associated with the possible dissociation channels, as shown in Table S2 (see Supplementary Material). Based on our calculations, the most favorable direct dissociation pathway for [${\mathrm{Ng}}_2{\mathrm{H]}}^{+}$ complexes involves the separation into a neutral Ng atom and ${\mathrm{NgH}}^{+}$ molecular cation. The relative energy associated with this process is $4633.10$ for He-complex and $5556.29$ for Ne-complex (see Figure 1). The next lowest dissociation energy corresponds to the loss of the proton, which requires an energy five times greater than predecessor dissociation channel. In this case, a competition is observed between the formation of ${\mathrm{Ng}}_2$ and the dissociation in 2Ng, with an energy difference of less than $1650$ in both scenarios. For [${\mathrm{He}}_2{\mathrm{H]}}^{+}$, the most favorable channel is the formation of ${\mathrm{Ng}}_2$, while in the case of [${\mathrm{Ne}}_2{\mathrm{H]}}^{+}$ the opposite mechanism is favored. A similar behavior is presented in the dissociation mechanism by loss of a hydrogen atom, with an energy twenty-three and fifteen times greater than the most favorable dissociation mechanism for the He-complex and Ne-complex, respectively. However, in this case, the energy difference between both pathways is much smaller, being approximately $15$ for [${\mathrm{He}}_2{\mathrm{H]}}^{+}$ and $12$ for [${\mathrm{Ne}}_2{\mathrm{H]}}^{+}$. These significant energy difference strongly indicates that the predominant formation way for these complexes is through the generation of the ${\mathrm{NgH}}^{+}$ and Ng atom. Consequently, our study of their interaction will focus on their dissociation into this lowest-energy channel.

In this context, Table S3 (see Supplementary Material) presents the stepwise formation energies as Ng atoms bind to the proton, revealing that the first Ng atom binds more strongly than the second in all studied systems. The total formation energy at T=0 K for the [${\mathrm{He}}_2{\mathrm{H]}}^{+}$ and [${\mathrm{Ne}}_2{\mathrm{H]}}^{+}$ systems is obtained by summing the stepwise formation energies, with the first binding step playing a dominant role in determining the overall behavior. For the ${\mathrm{He}}_2{\mathrm{H}}^{+}$, in general terms, an excellent agreement with the data available in the literature is observed. Compared with data in ref. [33], the discrepancies of almost 6 kcal/mol in the total formation energy and 1.5 kcal/mol in the second step are mainly due to the inclusion of zero-point energy corrections, which were not considered here, while the difference of approximately 1 kcal/mol with ref. [28] is due to the use of a different computational method (MRDCI). In the case of the ${\mathrm{Ne}}_2{\mathrm{H}}^{+}$, the observed discrepancies are again due to the differences in the computational methods used [35,39,40]. Comparing with the most recent study [35] a difference of about 1 kcal/mol is observed, suggesting that for noble gas hydrides the accuracy of the calculations becomes increasingly dependent on the size of the basis used as the Ng atom becomes heavier, as already observed in ref. [59,79].

2.1.2. Training and Testing Datasets

We have used the Jacobi’s coordinates to sample the configuration space for both [${\mathrm{He}}_2{\mathrm{H]}}^{+}$ and [${\mathrm{Ne}}_2{\mathrm{H]}}^{+}$ cations. The Jacobi (r, R, $\theta$) coordinates are defined as, r the vector along the Ng-H bond distance, R the vector along the distance between Ng from the center of mass of the H, and $\theta$ the angle between the (r, R) vectors. On the basis of our previous works [60,79], we have employed the single and double excitation coupled cluster with perturbative triples (CCSD(T)) method using the AV5Z and AV6Z basis sets, followed by extrapolation to the complete basis sets (CBS) to obtain the energies. We have utilized the two-point single inverse power function introduced by Schwartz [80] ${\mathrm{E}}_n^{corr}$ = ${\mathrm{E}}_{CBS}^{corr}$ + ${\displaystyle \frac{A}{n^3}}$, with n = 5 and 6, in order to compute the correlation energy at the CBS limit, and then used it in the ${\mathrm{E}}_{CBS}{\mathrm{=E}}_{HF/AV6Z}{\mathrm{+E}}_{CBS}^{corr}$ calculation.

The size of the training datasets was up to 8550 and 13680 CCSD(T)/CBS[56] ab initio points for the [${\mathrm{He}}_2{\mathrm{H]}}^{+}$ and [${\mathrm{Ne}}_2{\mathrm{H]}}^{+}$ molecular cations, respectively. The [${\mathrm{He}}_2{\mathrm{H]}}^{+}$ dataset was distributed in a grid of 10-points in the r coordinate ranging between 0.60 and 1.95 Å, 45-points in the R distance from 1.5 to 7.0 Å, and 19-points in the angular $\theta$ coordinate between 0 and ${180}^{\circ }$, while for the [${\mathrm{Ne}}_2{\mathrm{H]}}^{+}$ the grid is ranging from 0.75 to 2.25 Å, and 1.9 to 7.0 Å in r and R coordinates, respectively. For both molecules, we employed eight different sets of training data with in total 1000(1900)/1500(2850), 1800(3420)/2700(5130), 2500(4750)/4000(7600) and 4500(8550)/7200(13680) points by choosing 4 or 10/6 or 16, 25 or 45, and 10 or 19 grid points in the r, R and $\theta$ descriptors, respectively, for the [${\mathrm{He}}_2{\mathrm{H]}}^{+}$/[${\mathrm{Ne}}_2{\mathrm{H]}}^{+}$, keeping the same coordinate range in all of them. Once the generation of a high-quality dataset is completed, we then proceed with their organization, by splitting the whole dataset into training and test sets, utilized then in the training and evaluation processes, respectively. In Figure 2, we show the distributions of the cumulative number of configurations as a function of energy for the previously mentioned training sets, as well as, 11400 and 1840 additional testing dataset configurations of He- and Ne-complexes. We have also show distribution in energy of 3800 and 6080 randomly chosen data for each [${\mathrm{He}}_2{\mathrm{H]}}^{+}$ and [${\mathrm{Ne}}_2{\mathrm{H]}}^{+}$, respectively.

2.2. Potential Energy Surface Representations: Topology and Quality

2.2.1. RKHS ML-PESs Methodology

As mentioned, we choose the RKHS method proposed by Ho and Rabitz [68] to represent the PESs. Briefly, $V\left(x\right)$ represents the potential energy function, with $V_i$ denoting the known potential energies at specific molecular configurations $x_i$, then, by applying the representer theorem to a general functional relationship, $V\left(x\right)$ can be optimally approximated as a linear combination of appropriate functions: $\tilde{V}\left(x\right)\approx {\tilde{V}}_i={\sum }_{i=1}^{\mathcal{N}}{C}_{i}K(x,x_i)$, being $C_i$ coefficients and $K(x,x^{\prime })$ a kernel function. The coefficients $C_i$ are determined by solving the linear equation $\mathbf{C}={\mathbf{K}}^{-1}\mathbf{y}$, where $\mathbf{y}=[V_1,V_2,\dots ,V_{\mathcal{N}}]$ represents the input data, $\mathbf{C}=[C_1,\dots ,C_{\mathcal{N}}]$ is the vector of coefficients, and $\mathbf{K}$ is the $\mathcal{N}\times \mathcal{N}$ reproducing kernel matrix. Furthermore, the multidimensional kernel function $K(x,x^{\prime })$ can be expressed as a direct product: $K(x,x^{\prime })={\prod }_{d=1}^{D}k(x,x^{\prime })$, where D is the dimensionality of the kernel K, and k represents the one-dimensional (1D) kernels. Several 1D kernel functions for different purposes are available in the literature, and in this study, we have employed the $k_1^{n,m}$, $k_2^{n,m}$ and $k_3$ 1D reproducing kernel functions for the diatom distance-like (r), the Ng-diatom distance-like (R), and angle-like ($\theta$) variables, respectively. In the above equation, the reduced coordinate is $z=cos\theta$, whereas, $N_r$, $N_R$ and $N_{\theta }$ are the number of ab initio calculated points in each coordinate. The reproducing distance-like kernel functions $k_1$ and$k_2$ are represented by $k_{1,2}^{n,m}=n^2{\chi }_{>}^{-(m+1)}B{(m+1,n) }_2{F}_1(-n+1,m+1;n+m+1;{\displaystyle \frac{{\chi }_{<}}{{\chi }_{>}}})$, with $\chi$= r or R, respectively, ${\chi }_{>}$ and ${\chi }_{<}$ are the larger and smaller of the $\chi$, respectively, while the angle-like kernel function is given as $k_3(y,y^{\prime })={\sum }_l{\displaystyle \frac{(2ll+1)}{2}}{P}_l\left(y\right)P_l\left(y^{\prime }\right)$. The n and m superscripts refer to the order of smoothness of the function and its asymptotic behavior at large distances, with $n=2$, and m=3 accounting the $R^{-4}$ leading dispersion interaction between the Ng atom and ${\mathrm{NgH}}^{+}$ molecular ion. B is the beta function, ${}_{2}F_1$ is the Gauss hyper-geometric function, and ${\mathrm{P}}_{ll}$ the Legendre polynomials with $ll$=0-18. In turn, the potential form is given by

$$V(r,R,\theta )=\sum _{i=1}^{N_r}\sum _{j=1}^{N_R}\sum _{k=1}^{N_{\theta }}{C}_{ijk}{k}_1^{n,m}(r_i,r)k_2(R_j,R)k_3({\theta }_k,\theta ).$$

(1)

2.2.2. Validation of the RKHS ML-PES Models

A three-step validation protocol previously developed [60] was employed as follows: First, each PES’s model is validated on the training data to ensure that the results are consistent and reproducible with the information used during training. Then, potential values are generated at points outside the training range, using the test sets, to assess the ability of model to generalize to new and unseen data. Finally, regression diagnostics are monitored, here we will focus on the root-mean-square error (RMSE) and mean-absolute error (MAE), and it is verified that it does not exceed a predefined threshold. If the error exceeds this limit, more training data are incorporated to improve the model’s accuracy.

Figure 3 and Figure 4 illustrates the systematic improvement in the RMSE values of the RKHS ML-PES model as the number of training data increases. To select the best-performing RKHS PES for each complex, we employed a hold-out cross-validation scheme among the eight trained models. The RMSE values were computed considering the total number of testing data for [${\mathrm{Ng}}_2{\mathrm{H]}}^{+}$ cations. In the case of the He-complex, we averaged the RMSE over 3800 randomly selected points from the total 11400 test data, while for the Ne-complex, we used 6080 random points out of 18240. This allows us to assess the quality of the RKHS PES as a function of the training data size, with models trained on 4500/8550, 2500/4750, 1800/3420, and 1000/1900 or 7200/13680, 4000/7600, 2700/5130 and 1500/2850 points for both He or Ne-complexes (shown in the left and right plots of Figure S1 in Supplementary Material, respectively). As expected, RMSE values decrease as the dataset size increases. Additionally, we observed that RMSE values are highly sensitive to the sampling in the $\theta$ coordinate. Specifically, models trained on 4500/7200 (10,45,10/16,45,10), 2500/4000 (10,25,10/16,25,10), 1800/2700 (4,4,10/6,45,10), and 1000/1500 (4,25,10/6,25,10) exhibit significantly higher RMSE values compared to models trained on 8500/13680 (10,45,19/16,45,19), 4750/7600 (10,25,19/16,25,19), 3420/5130 (4,45,19/6,45,19), and 1900/2850 (4,25,19/6,25,19), where the RMSE values are significantly lower. The computed RMSE values for all RKHS ML-PES models, considering the total 11400/18240 test data, the 8866/17015 datasets with energies below the dissociation threshold, and the two corresponding randomly selected testing sets 3800(3713)/6080(5563), follow the same trend.

The correlation plots demonstrated the effectiveness of the selected RKHS ML-PES model by comparing its performance against 8500/13680 training and 11400/18240 testing data points for He/Ne molecular cations, covering both attractive and repulsive regions of the potential. The RMSE values, shown as a function of energy ranges (see upper panels of Figure 3 and Figure 4), indicate that outside the training region, they remain at 22/$23$, even for dissociation energies exceeding 8866/17015 configurations (see lower panels of Figure 3 and Figure 4). Additionally, the total mean absolute error (MAE) remains below 0.004% across 11400/18240 configurations.

When we analyze the overall quality of the RKHS ML-PES developed in this work, one can be seen that the amount of training data required for their construction is significantly lower compared to previous studies available in the literature for both complexes analyzed, such as those reported in [38] and [35]. For the case of [${\mathrm{He}}_2{\mathrm{H]}}^{+}$, there are two PESs in which RMSE values with respect to the ab initio data of 87.44 and $17.1$ are reported, references [34] and [38], respectively. The last one of these models used 15682 points with energies below $10000$, which contrasts with the non-error of the RKHS ML-PES presented here. In the case of [${\mathrm{Ne}}_2{\mathrm{H]}}^{+}$, the study by Koner et al. [35] reports an RMSE of $9.09$, based on more than 23000 potential energy values computed ab initio. In contrast with all of this, the RKHS ML-PESs developed in this work require a considerably smaller training set, demonstrating the efficiency of the kernel-based approach in the study cases. This approach significantly reduces the amount of data required to build a high-precision PES compared to traditional analytical methods, as evidenced by recent studies.

Apart from the typical quantitative error analysis of the RKHS ML-PES, we also proceed to analyze the behavior of the hydrides cations under study through their potential curves plots along representative coordinates.

Thus, Figure 5 and Figure 6 present a comparison of the RKHS ML-PES (in solid lines) and the CCSD(T)/CBS[56] energies (in circles symbols) as a function of r (upper panels), R (middle panels) and $\theta$ (lower panels). In the upper panels we show the effect on the interaction of the shared proton motion between the two Ng atoms for linear (upper left panels) and antilinear (upper right panels) configurations of the [${\mathrm{He}}_2{\mathrm{H]}}^{+}$ and [${\mathrm{Ne}}_2{\mathrm{H]}}^{+}$ complexes, respectively, keeping fixed the R distances at the range of 1.5 to 2.7 Å. In the middle panels the potential curves are plotted for fixed $\theta$ from ${\mathrm{0-180}}^{\circ }$ for r= 0.90/1.15 Å (left panels) and 1.35/1.55 Å (right panels) as a function of R, while minimum energy path (MEP) plots are shown in the lower panels as a function of $\theta$ optimizing both r and R coordinates. As it is seen the RKHS ML-PES smoothly reproduces the ab initio data for all analyzed cases. Additionally, the full 3D RKHS ML-potential for [${\mathrm{Ng}}_2{\mathrm{H]}}^{+}$ cations are also presented using two-dimensional contour plots in Figure S2 (see Supplementary material) in the ($\theta$, R) and (r, R) planes, over from -$24000$ to -$15000$ for both cases.

In this manner, one can observe that the global minimum corresponds to the linear [${\mathrm{Ng-H-Ng]}}^{+}$ configuration, while the local minimum coincides with the antilinear [${\mathrm{Ng-Ng-H]}}^{+}$ configuration in both cases. For the [${\mathrm{He}}_2{\mathrm{H]}}^{+}$ system, the well depth of the global minimum is equal to $21100.37$ in the RKHS ML-PES and is located at R=1.69 Å and r=0.95 Å, whereas the second minimum presents an energy of -$16778.80$ with R=2.27 Å. Regarding the [${\mathrm{Ne}}_2{\mathrm{H]}}^{+}$ complex, the well depth is equal to $24001.01$ in the ML-PES and is located at R=2.23 Å and r=1.15 Å, whereas the local minimum presents an energy of -$18728.80$ with R=2.69 Å. Moreover, the potential energy barrier between the global and local potential minima at R=2.27 Å with an energy of -$16759.10$ for [${\mathrm{He}}_2{\mathrm{H]}}^{+}$ and at R=2.70 Å with an energy of -$18724.95$ for [${\mathrm{Ne}}_2{\mathrm{H]}}^{+}$.

2.3. Bound-State Quantum Calculations: Vibrational Spectral Bands Assignment

We employed the MCTDH (multi-configuration time-dependent Hartree) method, as implemented in the MCTDH package of codes [81]. We solve the time-dependent Schrödinger equation [82,83], $i{\Psi }^{J,K}=H_{K,K}^J{\Psi }^{J,K}+H_{K,K\pm 1}^J{\Psi }^{J,K\pm 1}$, where ${\Psi }^J(r,R,\theta ,\varphi )={\displaystyle \frac{1}{\sqrt{2\pi }}}{\sum }_K{\Psi }^{JK}(r,R,\theta ){exp}^{iK\varphi }$, and the Hamiltonian terms being, $H_{K,K\pm 1}^J=-{\displaystyle \frac{\sqrt{J(J+1)-K(K\pm 1)}}{2{\mu }_1{R}^2}}$$\left(\mp {\displaystyle \frac{\partial }{\partial \theta }}-Kcot\theta \right)$ and $H_{K,K}^J=-{\displaystyle \frac{1}{2{\mu }_1}}{\displaystyle \frac{{\partial }^2}{\partial R^2}}-{\displaystyle \frac{1}{2{\mu }_2}}{\displaystyle \frac{{\partial }^2}{\partial r^2}}-{\displaystyle \frac{1}{2{\mu }_2{r}^{2}sin\theta }}\left({\displaystyle \frac{\partial }{\partial \theta }}sin\theta {\displaystyle \frac{\partial }{\partial \theta }}-{\displaystyle \frac{K^2}{sin\theta }}\right)+{\displaystyle \frac{1}{2{\mu }_1{R}^2}}\left(J(J+1)-2K^2-{\displaystyle \frac{1}{sin\theta }}\left({\displaystyle \frac{\partial }{\partial \theta }}sin\theta {\displaystyle \frac{\partial }{\partial \theta }}-{\displaystyle \frac{K^2}{sin\theta }}\right)\right)+V(r,R,\theta )$.

The ($r,R,\theta$) Jacobi coordinates are chosen with the R being parallel to the Z body-fixed (BF) axis and the r coordinate lying in the XZ-plane. The azimuthal Euler angle $\varphi$ specifying the orientation of the BF around R. The total angular momentum $J=L+j$ is conserved, with L and j being the orbital and diatomic angular momentum operators, respectively, and its projection along the Z, $J_Z$=K, being the same of that of $j_Z$, as $l_Z$=0 by construction. The ${\mu }_1={\displaystyle \frac{m_{Ng}(m_{Ng}+m{}_H)}{M}}$ and ${\mu }_2={\displaystyle \frac{m_{Ng}{m}_H}{m_{Ng}+m_H}}$ are the reduced masses, with $M=2m_{Ng}+m_H$ the total mass of ${\mathrm{Ng}}_2{\mathrm{H}}^{+}$ and $m_{Ng}$, $m_H$ the mass of Ng=He or Ne and H or D atoms (isotopes), respectively. The $V(r,R,\theta )$ term is the RKHS ML-PES for each ${\mathrm{Ng}}_2{\mathrm{H}}^{+}$ cation, and within the MCTDH method [84,85,86] the POTFIT algorithm [87] was employed to transform these terms as a sum of products of single-particle operators, namely natural potentials (NPs). The contraction is constructed over the angular coordinate $\theta$, and 25 NPs are used in the POTFIT calculations for each r and R coordinates. In these calculations we considered the relevant region below the ${\mathrm{Ng+Ng+H}}^{+}$ dissociation threshold with a root mean square error (RMSE) of the fit in this region smaller than $2.5$. The wavefunction is also expanded in a sum of products of time-dependent basis, $\psi (Q,t)$, namely single-particle functions (SPFs), as ${\Psi }^{JK}={\sum }_{nkl}{A}_{nkl}^J\left(t\right){\psi }_n^{\left(1\right)}(r,t){\psi }_k^{\left(2\right)}(R,t){\psi }_l^{\left(3\right)}(\theta ,K,t)$, with $A_{nkl}^J$ being the time-dependent coefficients and Q the nuclear coordinates. Both the wavefunction and Hamiltonian operator are represented on a grid of the ($r,R,\theta$) coordinates ranging in the configuration space of interest, with the variable K taking only integer values. We have run several convergence test, and we have found that convergence achieved using 51 harmonic oscillator (HO) DVR functions in each r and R coordinates in the range of [0.6,1.95]/[0.7,2.25] and [1.35,3.5]/[1.8,4.0] Å, respectively, while for the angular coordinate $\theta$ we used $55$ DVR functions in the [0,2$\pi$] interval.

The improved relaxation ($IR$) and the block improved relaxation ($BIR$) methods [88,89], as implemented in the MCTDH package [81], were used to calculate the ground and excited rovibrational states. We have used a set of (18,18,18) SPFs in each coordinate in the $IR$ and $BIR$ MCTDH calculations, and we calculated in total up to 20 and 50/70 rotational-vibrational states considering the H and D isotopes in the ${\mathrm{He}}_2{\mathrm{H}}^{+}{\mathrm{/D}}^{+}$ and ${\mathrm{Ne}}_2{\mathrm{H}}^{+}{\mathrm{/D}}^{+}$, respectively. All calculated J=0, 1 and 2 states are given in the Supplementary Material, while all vibrational energy levels and their assignment (when possible) are listed in Table 1 and Table 2 and also shown in Figure 5 and Figure 6 for the ${\mathrm{He}}_2{\mathrm{H}}^{+}$ and ${\mathrm{Ne}}_2{\mathrm{H}}^{+}$, respectively. Comparisons with values reported in previously reported theoretical studies [34,35,37] are included in Table 1 and Table 2.

The anharmonic ZPEs are 2310.5/1802.55 and 1933.1/$1384.65$ for the ${\mathrm{He}}_2{\mathrm{H}}^{+}{\mathrm{/Ne}}_2{\mathrm{H}}^{+}$ and ${\mathrm{He}}_2{\mathrm{D}}^{+}{\mathrm{/Ne}}_2{\mathrm{D}}^{+}$ , respectively, with deuteration reducing the ZPE values by 377.8/$417.9$. By comparing with the ZPEs previously reported, we found differences of 29 and $143$ for the ${\mathrm{He}}_2{\mathrm{H}}^{+}$, 85 and $172$ for the ${\mathrm{He}}_2{\mathrm{D}}^{+}$ [34,37], and $34$ [35] for the ${\mathrm{Ne}}_2{\mathrm{H}}^{+}$. One can also observe that the present ML-PESs predict 5/10 bound vibrationally excited levels for the ${\mathrm{He}}_2{\mathrm{H}}^{+}{\mathrm{/He}}_2{\mathrm{D}}^{+}$, with energies that are significantly different from previous results available for the ${\mathrm{He}}_2{\mathrm{H}}^{+}$ [34,37]. In turn, 34/70 bound vibrationally excited levels are obtained for the ${\mathrm{Ne}}_2{\mathrm{H}}^{+}{\mathrm{/Ne}}_2{\mathrm{D}}^{+}$, whose energies being in a better accord with the recently 22 reported ones [35] for the ${\mathrm{Ne}}_2{\mathrm{H}}^{+}$. Figure 7 and Table 2 display 2D probability distributions of various rovibrational ${\nu }_1,{\nu }_2^J,{\nu }_3$ states in the (r,R) or ($\theta$,R) or ($\theta$,r) planes, where ${\nu }_1$ is the anharmonic symmetric Ng–H⁺ stretch, ${\nu }_2$ the anharmonic doubly degenerate proton bend vibration, and ${\nu }_3$ the anharmonic asymmetric Ng–H⁺ stretching or shared-proton mode. The assignment of each energy level has been extracted by analyzing the corresponding nodal structure of the probability density functions (see corresponding plots in Figure 7 and Table 2). As shown, the ground vibrational state (0,0,0) has zero nodes, the (1,0,0) has one node in R, the (0,0,1) shows one node in r, while the (1,0,1) has one node in each r and R coordinate. In turn, the (0,$1^1$,0) and (0,2,0) density functions exhibit bending excitations with one and two nodes in $\theta$ coordinate, respectively.

Regarding the ${\nu }_1$ and ${\nu }_3$ bands, our results predict values of 962.1 and $1300.4$ for the ${\mathrm{He}}_2{\mathrm{H}}^{+}$, with ${\nu }_1$ closely matching previously reported theoretical values of 963.6 and $954.2$ [34,37], while the ${\nu }_2$ is found at $886.7$. A comparison with recent experimental data [8,22,23] for specific rovibrational bands of ${\mathrm{He}}_2{\mathrm{H}}^{+}$ is presented in Figure 7. Notably, the ${\nu }_3$ band has been observed at $1315.8444$, along with seven of its rovibrational lines [8], including the (0,0,$1^1$) transition at $1320.1882$. Additional estimates reported for the ${\nu }_2$ bending mode at $874.9$ or the ${\nu }_1$ at $956.6$ or the combination band ${\nu }_1+{\nu }_3$ at $2057.9$ [8,22,23], while for reference the ${\nu }_2$ and ${\nu }_3$ fundamentals have been estimated at 673 and $1014$ for the ${\mathrm{He}}_2{\mathrm{D}}^{+}$ [14]. Our predictions underestimate the ${\nu }_3$ and the ${\nu }_1+{\nu }_3$ band by 15.4 and $7.1$, respectively, while overestimating the ${\nu }_1$ and ${\nu }_2$ by 11.8 and $5.5$, respectively. In the case of the ${\mathrm{He}}_2{\mathrm{D}}^{+}$ our calculations predicted the three fundamental bands at 907.6, 717.9 and $1080.4$.

In turn, for the ${\mathrm{Ne}}_2{\mathrm{H}}^{+}$ the ${\nu }_1$ and ${\nu }_3$ vibrational bands are found at 475.5 and $1418.3$, with the ${\nu }_3$ band comparing well to the reported value of $1432$ [40]. Additionally, the ${\nu }_1+{\nu }_3$, 2${\nu }_1+{\nu }_3$ and 3${\nu }_1+{\nu }_3$ combination bands are predicted at 1807.3, 2163.8 and $2485.3$, respectively, in close agreement with the 1814, 2182 and $2541$ values obtained from a coupled anharmonic model [40]. To date no experimental measurements have been reported for the ${\mathrm{Ne}}_2{\mathrm{H}}^{+}{\mathrm{/Ne}}_2{\mathrm{D}}^{+}$ species, and thus the present predictions await experimental verification.

3. Conclusions

The present study focuses on the computational spectroscopic characterization of the proton-bound light noble-gas ${\mathrm{Ng}}_2{\mathrm{H}}^{+}$ cations. Our results were obtained from high-level and well-converged ab initio electronic structure calculations, and newly trained RKHS ML-PESs were developed using high-quality CCSD(T)/CSB[56] set of data; followed by nuclear quantum bound state computations within the MCTDH framework employing the ML-PESs. Rovibrational states were calculated and assignments for several vibrational bands were provided. Comparisons of the zero-point energies, fundamental vibrational frequencies, and combination bands with previously reported theoretical estimates and recent experimental measurements on ${\mathrm{He}}_2{\mathrm{H}}^{+}$ were also presented and discussed. We have found that vibrational fundamentals bands are in close agreement with the experimental values reported for the ${\mathrm{He}}_2{\mathrm{H}}^{+}$, while our results for the ${\mathrm{Ne}}_2{\mathrm{H}}^{+}$ provide detailed and precise spectroscopic predictions that may serve as benchmarks for future laboratory or astrophysical detections.

Looking ahead, the quantum methodologies and insights developed in this work can be extended to the heavier proton-bound noble gas complexes, broadening the exploration of spectral features of this intriguing class of molecular ions. The success of the ML-PES strategy also opens avenues for modeling larger and more complex weakly bound systems, bridging the gap between computational predictions and experimental efforts - an ongoing challenge for the computational spectroscopy. We anticipate that the present results will stimulate further experimental studies, and will contribute to a deeper understanding of such noble gas cations in the ISM.