## Abstract

The high-lying vibrational states of the magnesium dimer (Mg_{2}), which has been recognized as an important system in studies of ultracold and collisional phenomena, have eluded experimental characterization for half a century. Until now, only the first 14 vibrational states of Mg_{2} have been experimentally resolved, although it has been suggested that the ground-state potential may support five additional levels. Here, we present highly accurate ab initio potential energy curves based on state-of-the-art coupled-cluster and full configuration interaction computations for the ground and excited electronic states involved in the experimental investigations of Mg_{2}. Our ground-state potential unambiguously confirms the existence of 19 vibrational levels, with ~1 cm^{−1} root mean square deviation between the calculated rovibrational term values and the available experimental and experimentally derived data. Our computations reproduce the latest laser-induced fluorescence spectrum and provide guidance for the experimental detection of the previously unresolved vibrational levels.

## INTRODUCTION

The weakly bound alkaline-earth dimers (AE_{2}) have emerged as probes of fundamental physics relevant to ultracold collisions (*1*), doped helium nanodroplets (*2*), coherent control of binary reactions (*3*), and even fields rarely associated with molecular science, such as optical lattice clocks (*4*) and quantum gravity (*5*). The magnesium dimer is especially important, since it has several desirable characteristics that can be useful in the above applications, such as the absence of hyperfine structure in the most abundant ^{24}Mg isotope that facilitates the analysis of binary collisions involving laser-cooled and trapped atoms, it helps us understand heavier AE_{2} diatomics, and, unlike its lighter Be_{2} analog, it is nontoxic (*6*). Unfortunately, the status of Mg_{2} as a prototype heavier AE_{2} species is complicated by the fact that its high-lying vibrational levels and, consequently, the long-range part of its ground-state potential energy curve (PEC) have eluded experimental characterization for half a century. In this regard, the magnesium dimer is even more challenging than its celebrated beryllium counterpart, whose elusive 12th vibrational level near the dissociation threshold (*7*, *8*), which we also found in (*9*), was confirmed in 2014 (*10*) after reanalyzing the spectra obtained in stimulated emission pumping experiments (*11*).

Experimentally, probing vibrational manifold of the magnesium dimer in its ground, _{2}, being a homonuclear diatomic, is infrared inactive. The first high-resolution photoabsorption spectra of Mg_{2}, corresponding to a transition from the ground state to the electronically excited *12*). Their spectroscopic analysis resulted in 285 *G*(*v*″, *J*″) and 656 *G*(*v*′, *J*′) rovibrational term values of ^{24}Mg_{2} involving 13 (*v*″ = 0 to 12) *v*′ = 1 to 24) *v*, and rotational, *J*, quantum numbers in the ground electronic state are designated by a double prime, whereas those corresponding to the excited *13*–*16*) *v*″ = 12 level about 25 cm^{−1} below the dissociation threshold, pointing to the existence of extra vibrational states with *v*″ > 12. It did not take long to detect one of such states. In 1973, Li and Stwalley (*17*) identified *v*″ = 13 level in the spectra reported in (*12*). They accomplished this by extending the original RKR PEC of Balfour and Douglas to the asymptotic region beyond 7.16 Å using theoretical values of *C*_{6} and *C*_{8} van der Waals coefficients (*18*, *19*). The resulting PEC supported 19 vibrational levels, i.e., five levels more than what was observed experimentally (*17*). Four decades later, in an effort to characterize states with *v*″ > 13, Knöckel *et al.* (*20*, *21*) examined the *22*). They improved and expanded the original ^{24}Mg_{2} dataset of Balfour and Douglas by reporting a total of 333 *G*(*v*″, *J*″) and 1,351 *G*(*v*′, *J*′) rovibrational term values involving *v*″ = 0 to 13 and *v*′ = 1 to 46, respectively, and constructed a few experimentally derived analytical forms of the *C*_{6} (*23*), *C*_{8} (*24*), and *C*_{10} (*24*) coefficients, which support the discrete spectral data in the 3.27 to 8.33 Å range (*20*). Although these refined PECs supported 19 ^{24}Mg_{2} vibrational levels, reinforcing the initial prediction of Li and Stwalley (*17*), Knöckel *et al.* (*20*) were unable to identify *v*″ > 13 in their LIF spectra.

Typically, high-lying vibrational states near dissociation constitute a small fraction of the entire vibrational manifold, but this is not the case for the weakly bound magnesium dimer, which has a shallow minimum on the ground-state PEC at *r*_{e} = 3.89039 Å (*20*) and a tiny dissociation energy *D*_{e} of 430.472(500) cm^{−1} (*20*, *21*). If the five extra levels, which have been speculated about, truly existed, they would represent more than a quarter of the entire vibrational manifold in the ground electronic state. Furthermore, without precise knowledge of the ground-state PEC of Mg_{2}, especially its long-range part that determines the positions of the high-lying vibrational states near the dissociation threshold, one cannot accurately interpret the aforementioned ultracold and collisional phenomena involving interacting magnesium atoms. It is intriguing why a seemingly docile main group diatomic continues to challenge state-of-the-art spectroscopic techniques. The experimental difficulties in detecting the elusive *v*″ > 13 states of the magnesium dimer originate from several factors, including small energy gaps between high-lying vibrations that are comparable to rotational spacings (*12*, *25*), resulting in overlapping spectral lines, and unfavorable signal-to-noise ratio in the existing LIF spectra (*20*). Rotational effects complicate the situation even more, since, in addition to affecting line intensities (*20*, *22*, *25*), they may render the high-lying vibrational states of Mg_{2} unbound. All of these and similar difficulties prompted Knöckel *et al.* (*20*, *21*) to conclude that experimental work alone is insufficient and that accurate theoretical calculations are needed to guide further analysis of the ground-state PEC and rovibrational states of Mg_{2}, especially the elusive *v*″ > 13 levels near the dissociation threshold.

Unfortunately, there have only been a handful of theoretical investigations attempting to determine the entire vibrational manifold of the magnesium dimer. This is, at least in part, related to the intrinsic complexity of the underlying electronic structure and difficulties with obtaining an accurate representation of the ground-state PEC using purely ab initio quantum-chemical means. At the Hartree-Fock theory level, which neglects electron correlation and dispersion interactions, Mg_{2} remains unbound. As demonstrated in this work, one needs to go to much higher theory levels, incorporate high-order many-electron correlation effects, including valence as well as inner-shell electrons, and use large, carefully calibrated, one-electron basis sets to accurately capture the relevant physics and obtain a reliable description of the *26*) for a detailed discussion and historical account, including references to the earlier quantum chemistry computations for the magnesium dimer]. Ab initio quantum mechanical calculations for the

The initial theoretical estimates of the number of vibrational states supported by the *27*), while the more recent ab initio quantum chemistry computations based on the various levels of coupled-cluster (CC) theory (*28*), reported in (*26*, *29*), suggested that the highest vibrational level of ^{24}Mg_{2} is *v*″ = 18. Among the previous theoretical studies, only Amaran *et al.* (*29*) considered the *12*). Furthermore, as demonstrated in our recent benchmark study (*26*), where a large number of CC methods were tested using the ^{24}Mg_{2} as examples, and consistent with the earlier calculations (*30*, *31*), the popular CCSD(T) approximation (*32*) exploited in (*29*) could not possibly produce the small RMSD value reported in (*29*), of 1.3 cm^{−1}, for the rovibrational manifold of Mg_{2} in its ground electronic state; the value on the order of a dozen cm^{−1} would be more appropriate (*26*) (see fig. S1). Similar remarks apply to the *29*) using the low-level variant of the linear-response CC theory (*33*), resulting in noticeable deviations from the experimentally derived *21*). To simulate and properly interpret the *20*) using purely theoretical means, one needs much higher accuracy levels in the computations of line positions and robust information about line intensities, which has not been obtained in the previous quantum chemistry studies.

The call for a reliable ab initio computation of the ground-state PEC and rovibrational states of Mg_{2}, including the *v*″ > 13 levels that have eluded experimentalists for decades, expressed by Knöckel *et al.* (*20*, *21*), is answered in the present work. We report the highly accurate PECs for the ground, _{2} (*12*, *20*–*22*), obtained with state-of-the-art ab initio quantum chemistry, and use them to determine the corresponding rovibrational manifolds. Consistent with the conclusions of our recent benchmark study (*26*), to obtain a highly accurate representation of the ground-state PEC, we combine the numerically exact description of the valence electron correlation effects, provided by the full configuration interaction (CI) approach, with the nearly exact description of subvalence correlations involving all electrons but the 1s shells of Mg atoms offered by the CC theory with a full treatment of singly, doubly, and triply excited clusters, abbreviated as CCSDT (*34*, *35*). Our computational protocol for the *26*), is similar, except that to capture subvalence electron correlation effects in this state we adopt one of the carefully chosen approximations to the equation-of-motion (EOM) CC theory (*36*) with singles, doubles, and triples (EOMCCSDT) (*37*, *38*) belonging to the completely renormalized CR-EOMCCSD(T) family (*39*), which is considerably more affordable than EOMCCSDT without substantial loss of accuracy. As in the case of the *12*, *17*, *20*, *21*), and to accurately describe the relevant many-electron correlation effects, the electronic structure calculations reported in this work rely on the carefully calibrated augmented polarized valence and weighted core-valence correlation-consistent bases of quadruple-ζ quality developed in (*40*), designated as aug-cc-pV(Q + *d*)Z and aug-cc-pwCVQZ, respectively. We used these basis sets in our earlier benchmark calculations for the ground state of the magnesium dimer (*26*) but only for the methods up to CCSDT, i.e., the important post-CCSDT electron correlation effects were treated in (*26*) with smaller, less saturated, bases.

To make our comparisons with the experiment more complete, for each of the two electronic potentials considered in this study, we examine both the most abundant ^{24}Mg_{2} species and the ^{24}Mg^{25}Mg, ^{24}Mg^{26}Mg, ^{25}Mg_{2}, ^{25}Mg^{26}Mg, and ^{26}Mg_{2} isotopologs (to our knowledge, rovibrational levels of the Mg_{2} species other than ^{24}Mg_{2} have not been calculated using ab initio potentials before). We combine the above information with the *20*, *21*), including line positions and the corresponding line intensities, as defined via the Einstein coefficients, and provide the long-awaited theoretical guidance for the possible experimental detection of the *v*', *J*') → *v″*, *J″*) rovibronic transitions involving the *v″* > 13 levels.

## RESULTS

The most essential numerical information, generated in the present study using the computational protocol described in Materials and Methods, is summarized in Figs. 1 to 3 and Tables 1 to 3. All of the numerical data supporting the content and conclusions of this work are included in the main text and compiled in the Supplementary Materials document and data files S1 and S2 attached to it. In describing and discussing our results, we begin with the PECs and rovibrational term values characterizing the *12*, *20*, *21*). Next, we compare the experimental LIF spectra reported in (*20*, *21*) with those resulting from our computations and suggest potential avenues for detection of the elusive *v*″ > 13 levels of the magnesium dimer. Auxiliary information, which complements the discussion in this section, including further comments on the accuracy and convergence characteristics of the computational protocol used in the present study, the effect of isotopic substitution on the calculated rovibrational term values, the discussion of the validity of the Franck-Condon analysis adopted in (*20*) to examine the LIF spectra reported in (*20*, *21*), and the lifetimes for predissociation by rotation characterizing quasi-bound rovibrational states supported by the

### PECs and rovibrational states

As shown in Table 1, our ab initio *D*_{e} and equilibrium bond length *r*_{e} of Mg_{2} (*20*, *21*) to within 0.9 cm^{−1} (0.2%) and 0.003 Å (0.07%), respectively. These high accuracies in describing *D*_{e} and *r*_{e} are reflected in our calculated rovibrational term values of ^{24}Mg_{2} and its isotopologs, which are in very good agreement with the available experimental information (*12*, *20*, *21*). As shown in the spreadsheets included in data file S1, the RMSDs characterizing our ab initio *G*(*v*″, *J*″) values for ^{24}Mg_{2} relative to their experimentally determined counterparts, reported in (*12*) for *v*″ < 13 and (*20*, *21*) for *v*″ < 14, are 1.1 cm^{−1}, when the spectroscopic data from (*12*) are used, and 1.5 cm^{−1}, when we rely on (*20*, *21*) instead. At the same time, the maximum unsigned errors in our calculated *G*(*v*″, *J*″) values relative to the experiment do not exceed ~2 cm^{−1}, even when the quasi-bound states above the potential asymptote arising from centrifugal barriers are considered. Although the experimental information about the *G*(*v*″, *J*″) values characterizing other Mg_{2} isotopologs is limited to ^{24}Mg^{25}Mg, ^{24}Mg^{26}Mg, and ^{26}Mg_{2} and includes very few *v*″ values (*20*, *21*), the RMSDs relative to the experiment resulting from our calculations are similarly small (1.0 cm^{−1} for ^{24}Mg^{25}Mg, 1.2 cm^{−1} for ^{24}Mg^{26}Mg, and 0.6 cm^{−1} for ^{26}Mg_{2}; cf. the Supplementary Materials).

Further insights into the quality of our ab initio calculations for the ground-state PEC can be obtained by comparing the resulting rovibrational term values with their counterparts determined using the most accurate, experimentally derived, analytical forms of the *20*). In the discussion below, we focus on the so-called X-representation of the ground-state PEC developed in (*20*), which the authors of (*20*) regard as a reference potential in their analyses (see Table 2). We recall that the X-representation of the ground-state PEC of the magnesium dimer was obtained by simultaneously fitting the *C*_{6} (*23*), *C*_{8} (*24*), and *C*_{10} (*24*) coefficients. As shown in Table 2, our ab initio *G*(*v*″, *J*″) energies characterizing the most abundant ^{24}Mg_{2} isotopolog are in very good agreement with those generated using the X-representation of the ground-state PEC developed in (*20*). When all of the rovibrational bound states supported by both potentials are considered, the RMSD and the maximum unsigned error characterizing our ab initio *G*(*v*″, *J*″) values for ^{24}Mg_{2} relative to their counterparts arising from the X-representation are 1.3 and 2.0 cm^{−1}, respectively. What is especially important in the context of the present study is that our ab initio ground-state PEC and the state-of-the-art analytical fit to the experimental data defining the X-representation, constructed in (*20*), bind the *v*″ = 18 level if the rotational quantum number *J*″ is not too high (see the discussion below).

The high quality of our calculated *G*(*v*″, *J*″) values and spacings between them, which can also be seen in Tables 1 and 2 and Fig. 1, allows us to comment on the existence of the *v*″ > 13 levels that have escaped experimental detection for decades. As already alluded to above and as shown in Table 2 and Fig. 1, our ab initio *20*), which for the most abundant ^{24}Mg_{2} isotopolog is 19 (see the Supplementary Materials for the information about the remaining Mg_{2} species). Table 1, which compares the rovibrational term values of ^{24}Mg_{2} resulting from our ab initio calculations for the representative rotational quantum numbers ranging from 0 to 80 with the available experimental data, shows that the elusive high-lying states with *v*″ > 13 quickly become unbound as *J*″ increases, so by the time *J*″ = 20, the *v*″ = 15 to 18 levels are no longer bound (see fig. S3 for a graphical representation of the *J*″ = 20, 40, 60, and 80 effective potentials including centrifugal barriers characterizing the rotating ^{24}Mg_{2} molecule, along with the corresponding vibrational wave functions and information about the lifetimes for predissociation by rotation associated with tunneling through centrifugal barriers characterizing quasi-bound states). In fact, according to our ab initio data compiled in the Supplementary Materials, the maximum rotational quantum number that allows for at least one bound rovibrational state decreases with *v*″, from *J*″ = 68 for *v*″ = 0 to *J*″ = 4 for *v*″ = 18, with all states becoming quasi-bound or unbound when *J*″ ≥ 70, when the most abundant ^{24}Mg_{2} isotopolog is considered. In general, as shown in fig. S3 and the lifetime data compiled in the Supplementary Materials, the mean lifetimes for predissociation by rotation characterizing quasi-bound states with a given *J*″ rapidly decrease as *v*″ becomes larger. They decrease equally fast when *J*″ increases and *v*″ is fixed. These observations imply that the spectroscopic detection of the high-lying vibrational states of Mg_{2} can only be achieved if the molecule does not rotate too fast (cf. Table 1 and fig. S3). We could not find any information regarding the timescales involved in the LIF experiments carried out by Knöckel *et al.* (*20*). However, a comparison of our ab initio–determined quasi-bound rovibrational states, including their energies and lifetimes compiled in the Supplementary Materials, with the observed rovibronic transitions reported in the supplementary material of (*21*) suggests that the mean lifetimes for predissociation by rotation characterizing quasi-bound states seen in the experimentally resolved LIF spectral lines are on the order of 0.1 ns or longer.

As shown in Fig. 1, where we plot the wave functions of the high-lying, purely vibrational, states of ^{24}Mg_{2}, starting with the last experimentally observed *v*″ = 13 level, along with the *v*″ = 18 state, located only 0.2 cm^{−1} below the potential asymptote, is barely bound (see also Table 1). This makes the existence of an additional, *v*″ = 19, level for the most abundant isotopolog of the magnesium dimer unlikely. Further insights into the number of purely vibrational bound states of ^{24}Mg_{2} supported by the *G*(*v*″ + 1) − *G*(*v*″) energy differences, resulting from the ab initio calculations reported in this work and the experiment, as a function of *v*″ + ½ (the Birge-Sponer plot). Fitting the experimental data to a line, i.e., assuming a Morse potential, results in *v*″ = 16 being the last bound vibrational level of ^{24}Mg_{2}. Although the deviation from the Morse potential, as predicted by our ab initio calculations, is not as severe as in the case of Be_{2} (*11*), it is large enough to result in the *v*″ = 17 and 18 states becoming bound, emphasizing the importance of properly describing the long-range part of the PEC.

As shown in Table 1 and Fig. 1, the *G*(*v*″ + 1) − *G*(*v*″) vibrational spacings rapidly decrease with increasing *v*″, from 47.7 cm^{−1} or 68.6 K for *v*″ = 0 to 11.7 cm^{−1} or 16.8 K for *v*″ = 12, and to 0.8 cm^{−1} or 1.2 K for *v*″ = 17, when ^{24}Mg_{2} is considered. This means that at regular temperatures all vibrational levels of the magnesium dimer, which is a very weakly bound system, are substantially populated, making selective probing of the closely spaced higher-energy states, including those with *v*″ > 13, virtually impossible, since practically every molecular collision (e.g., with another dimer) may result in a superposition of many rovibrational states, with some breaking the dimer apart. At room temperature, for example, the cumulative population of the *v*″ > 13 states of ^{24}Mg_{2}, determined using the normalized Boltzmann distribution involving all rotationless levels bound by the *v*″ = 0, 13% for *v*″ = 1, and 10% for *v*″ = 2). The situation changes in the cold/ultracold regime, where the available thermal energies, which are on the order of millikelvin or even microkelvin, are much smaller than the vibrational spacings, even when the high-lying states with *v*″ > 13 near the dissociation threshold are considered, suppressing collisional effects and allowing one to probe the long-range part of the ground-state PEC, where the *v*″ > 13 states largely localize (cf. Fig. 1). This makes the accurate characterization of the *v*″ > 13 bound and quasi-bound states provided by the high-level ab initio calculations reported in this work relevant to the applications involving cold/ultracold Mg atoms separated by larger distances in magneto-optical traps [see, e.g., (*6*)].

The accuracy of our ab initio description of the more strongly bound *D*_{e} = 9414 cm^{−1} and *r*_{e} = 3.0825 Å (*21*); cf. Fig. 2 for the corresponding PEC], which we need to consider to simulate the LIF spectra, is consistent with that obtained for the weakly bound ground state. For example, the errors relative to the experiment (*21*) resulting from our calculations of the dissociation energy *D*_{e} and equilibrium bond length *r*_{e} are 0.91% (86 cm^{−1}) and 0.2% (0.006 Å), respectively (see the Supplementary Materials). This high accuracy of our ab initio ^{24}Mg_{2} *G*(*v*′, *J*′) values obtained in this work and their experimentally derived counterparts reported in (*12*, *21*). In particular, the RMSDs characterizing our rovibrational term values in the *12*) and Knöckel *et al.* (*21*) are only 3.2 and 4.5 cm^{−1}, respectively, which is a major improvement over the RMSD of 30 cm^{−1} reported in (*29*). They are similarly small for the rovibrational states supported by the ^{24}Mg^{25}Mg, ^{24}Mg^{26}Mg, and ^{26}Mg_{2} isotopologs examined in (*21*) (3.7, 4.1, and 4.1 cm^{−1}, respectively; cf. the Supplementary Materials). According to our ab initio calculations using the computational protocol described in Materials and Methods, the total number of vibrational states supported by the ^{24}Mg_{2} species is 169 (see data file S1).

### LIF: Ab initio theory versus experiment

The most compelling evidence for the predictive power of our ab initio electronic structure and rovibrational calculations is the nearly perfect reproduction of the experimental *20*, *21*), shown in Fig. 3 and Table 3, with further information provided in the Supplementary Materials. Figure 2 uses our calculated *20*), which is reproduced in Fig. 3A. This particular spectrum represents the fluorescence progression from the *v*′ = 3, *J*′ = 11) state of ^{24}Mg_{2}, populated by laser excitation from the *v*″ = 5, *J*″ = 10) state, to all accessible *v*″, *J*″) rovibrational levels, resulting in the P12/R10 doublets that correspond to *J*″ = 12 for the P branch and *J*″ = 10 for the R branch. Figure 3 and Table 3 compare the experimentally observed *v*′ = 3, *J*′ = 11) → *v*″, *J*″ = 10,12) transitions with the corresponding line positions (Fig. 3 and Table 3) and intensities (Fig. 3) resulting from our ab initio calculations. The only adjustment that we made to produce the theoretical LIF spectrum shown in Fig. 3 and Table 3 was a uniform shift of the entire *T*_{e} of 26,068.9 cm^{−1} (*21*) (see Materials and Methods for the details). Other than that, the theoretical LIF spectrum in Fig. 3 and Table 3 relies on the raw ab initio electronic structure and rovibrational data.

The notable agreement between the theoretical and experimental LIF spectra shown in Fig. 3A and Table 3, with differences in line positions not exceeding 1 to 1.5 cm^{−1} and with virtually identical intensity patterns, suggests that our predicted transition frequencies involving the elusive *v*″ > 13 states are very accurate, allowing us to provide guidance for their potential experimental detection in the future. Before discussing our suggestions in this regard, we note that owing to our ab initio calculations, we can now locate the previously unidentified P12/R10 doublets involving the *v*″ > 13 states within the experimental LIF spectrum reported in figure 3 of (*20*). Indeed, as shown in Fig. 3 and Table 3, the LIF spectrum corresponding to the *v*′ = 3, *J*′ = 11) → *v*″, *J*″ = 10,12) transitions contains the P12/R10 doublets involving the *v*″ = 0 to 16 states and the R10 line involving the *v*″ = 17 state. The *v*′ = 3, *J*′ = 11) → *v*″ = 17, *J*″ = 12) and *v*′ = 3, *J*′ = 11) → *v*″ = 18, *J*″ = 10,12) transitions are absent, since the *v*″ = 17, *J*″ = 12 and *v*″ = 18, *J*″ = 10 and 12 states are unbound (see the Supplementary Materials), but they could potentially be observed if one used different initial *v*′, *J*′) states (see the discussion below).

As one can see by inspecting data file S2 and Fig. 3, and consistent with the remarks made by Knöckel *et al.* in (*20*), the experimental detection of the P12/R10 doublets involving *v*″ > 13, when transitioning from the *v*′ = 3, *J′* = 11) state, was hindered by the unfavorable signal-to-noise ratio (transitions to the *v*″ = 16 and 17 states exhibit low Einstein coefficients) and the presence of overlapping lines outside the P12/R10 progression, originating from collisional relaxation effects (*20*) and having similar (*v*″ = 15) or higher (*v*″ = 14) intensities. To fully appreciate this, in Fig. 3B, we magnified the region of the LIF spectrum recorded in (*20*) that contains the calculated *v*′ = 3, *J*′ = 11) → *v*″ = 13 to 16, *J*″ = 10,12) and *v*′ = 3, *J*′ = 11) → *v*″ = 17, *J*″ = 10) transitions. As shown in Fig. 3 and Table 3, the identification of the P12/R10 doublets corresponding to the *v*′ = 3, *J*′ = 11) → *v*″ = 0 to 13, *J*″ = 10,12) transitions is unambiguous. The observed and calculated line positions and intensities and line intensity ratios within every doublet match each other very closely. Figure 3B demonstrates that the identification of the remaining doublets in the P12/R10 progression is much harder. On the basis of our ab initio work and taking into account the fact that our calculated line positions may be off by about 1 cm^{−1} (cf. Table 3), the *v*″ = 14 P12/R10 doublet, marked in Fig. 3B by the blue arrows originating from the *v*″ = 14 label, is largely hidden behind the higher-intensity feature that does not belong to the P12/R10 progression and that most likely originates from collisional relaxation (*20*). Because of our calculations, we can also point to the most likely location of the *v*″ = 15 P12/R10 doublet in the LIF spectrum recorded in (*20*) (see the blue arrows originating from the *v*″ = 15 label in Fig. 3B). Doing this without backing from the theory is virtually impossible due to the presence of other lines near the *v*′ = 3, *J*′ = 11) → *v*″ = 15, *J*″ = 10,12) transitions having similar intensities. As shown in Fig. 3B, the situation with the remaining *v*′ = 3, *J*′ = 11) → *v*″ = 16, *J*″ = 10,12) and *v*′ = 3, *J*′ = 11) → *v*″ = 17, *J*″ = 10) transitions is even worse, since they have very low Einstein coefficients that hide them in the noise.

Knöckel *et al.* (*20*) also suggested that the difficulties with detecting the P12/R10 doublets involving the *v*″ = 14 and 15 states, which have higher Franck-Condon factors than those characterizing the experimentally observed *v*′ = 3, *J*′ = 11) → *v*″ = 0, *J*″ = 10,12) transitions, might be related to the variation of the *r* and limitations of the Franck-Condon principle, but our ab initio calculations do not confirm this. As shown in the Supplementary Materials (see fig. S2), in the region of *r* values where the respective rovibrational wave functions, *v*″ = 14,15 and *J*″ = 10,12 for the *v*′ = 3 and *J*′ = 11 for the *r* = 2.2 to 100.0 Å region examined in this work is considered. In agreement with the analysis presented in (*20*), our calculated Franck-Condon factors for the *v*″ = 14 and 15 P12/R10 doublets are higher than those characterizing the analogous *v*″ = 0 lines, but, as demonstrated in data file S2, the same holds for the respective Einstein coefficients. Thus, it is the presence of densely spaced and overlapping lines outside the P12/R10 progression having similar or higher intensities than the *v*′ = 3, *J*′ = 11) → *v*″ = 14,15, *J*″ = 10,12) transitions that makes the experimental identification of the *v*″ = 14 and 15 P12/R10 doublets very hard.

### Theory-inspired avenues for detection of elusive states

In general, our ab initio calculations carried out in this work indicate that under the constraints of the LIF experiments reported in (*20*, *21*), where the authors populated the *v*′, *J*′) states with *v*′ = 1 to 46, the *v*″, *J*″) states with *v*″ = 14 to 18 cannot be realistically detected because of very small Franck-Condon factors and Einstein coefficients characterizing the corresponding *v*′, *J*′) → *v*″, *J*″) transitions (see data file S2). As shown in Fig. 1, the *v*″ = 14 to 18 states are predominantly localized in the long-range *r* = 8 to 16 Å region. At the same time, as illustrated in Fig. 2, the potential well characterizing the electronically excited *v*″, *J*″) states with *v*″ = 14 to 18 via fluorescence from *v*′, *J*′) levels with *v*′ ≫ 46.

In an effort to assist the experimental community in detecting the elusive *v*″ = 14 to 18 vibrational levels, we searched for the *v*′, *J*′) → *v*″ = 14 to 18, *J*″ = *J*′ ± 1) transitions in the most abundant isotopolog of the magnesium dimer, ^{24}Mg_{2}, that would result in spectral lines of maximum intensity based on the Einstein coefficients compiled in data file S2. To ensure the occurrence of allowed transitions involving the last, *v*″ = 18, level, which for ^{24}Mg_{2} becomes unbound when *J*″ > 4, we focused on the *J*″ values not exceeding 4, i.e., the fluorescence from the *v*′, *J*′) states with *J*′ = 1, 3, and 5. According to our calculations, the optimum *v*′ values for observing the *v*″ = 14 to 18, *J*″ ≤ 4 states via the LIF spectroscopy are in the neighborhood of *v*′ = 60, 66 to 69, and 74 to 84 for *v*″ = 14; 72 to 75 and 80 to 91 for *v*″ = 15; 79 to 82 and 88 to 100 for *v*″ = 16; 88, 89, and 97 to 111 for *v*″ = 17; and 109 to 129 for *v*″ = 18 [see data file S2 for the details of all allowed rovibronic transitions in ^{24}Mg_{2} involving the *v*″, *J*″ ≤ 4) → *v*′, *J*′) pump and *v*′, *J*′ = 1,3,5) → *v*″ = 14 to 18, *J*″ ≤ 4) fluorescence processes]. In determining these optimum *v*′ values, we chose the cutoff value of 1.0 × 10^{7} Hz in the Einstein coefficients, which is similar to the Einstein coefficients calculated for the most intense *v*″ = 5 P12/R10 doublet in the experimental LIF spectrum shown in figure 3 of (*20*), reproduced in Fig. 3A. Our predicted *v*′, *J*′ = 1,3,5) → *v*″ = 14 to 18, *J*″ ≤ 4) fluorescence frequencies resulting from the aforementioned optimum *v*′ ranges, which might allow one to detect the *v*″ = 14 to 18 states of ^{24}Mg_{2} via a suitably designed LIF experiment, are estimated at about 33,360, 33,740 to 33,910, and 34,150 to 34,530 cm^{−1} for *v*″ = 14; 34,050 to 34,190 and 34,390 to 34,710 cm^{−1} for *v*″ = 15; 34,350 to 34,460 and 34,640 to 34,880 cm^{−1} for *v*″ = 16; 34,630 to 34,660 and 34,830 to 35,000 cm^{−1} for *v*″ = 17; and 34,990 to 35,100 cm^{−1} for *v*″ = 18 [given the 86 cm^{−1} error in the calculated *D*_{e} characterizing the ^{−1} in our ^{24}Mg_{2} *G*(*v*′, *J*′) values relative to the spectroscopic data of (*12*, *21*), the above frequency ranges may have to be shifted by a dozen or so cm^{−1}].

## DISCUSSION

We used state-of-the-art ab initio quantum-mechanical methodologies to address a half-century-old enigma regarding the *v*″ = 14 to 18 vibrational states of the magnesium dimer. We provided the highly accurate ground-state PEC and rovibrational term values of ^{24}Mg_{2} and its less abundant ^{24}Mg^{25}Mg, ^{24}Mg^{26}Mg, ^{25}Mg_{2}, ^{25}Mg^{26}Mg, and ^{26}Mg_{2} isotopologs. We demonstrated that the ^{24}Mg_{2} up to *v*″ = 18, although the elusive *v*″ > 13 states become unbound as the rotational quantum number *J*″ increases, which contributes to difficulties with their experimental detection. We also obtained an accurate representation of the ^{24}Mg_{2}, and, with the help of the ab initio electronic transition dipole moment function, determined in this study as well, accurately simulated the LIF spectra recorded in (*20*, *21*), including line positions and intensities. Our work provides the long-awaited guidance for possible experimental identification of rovibronic transitions involving the *v*″ > 13 levels that have eluded scientists for five decades.

We hope that this study will fuel new spectroscopic investigations of the challenging Mg_{2} species and its heavier group 2 analogs, which are important in a variety of phenomena at the intersection of chemistry and atomic, molecular, and optical physics. A few years ago, ab initio calculations (*8*) combined with spectroscopic analyses (*7*, *10*) led to the discovery of the elusive 12th vibrational level of the beryllium dimer. By dealing with five similarly challenging states in a system three times larger than Be_{2}, we demonstrated that the predictive power of modern ab initio quantum chemistry is no longer limited to small few-electron species.

## MATERIALS AND METHODS

### Ab initio electronic structure calculations

The goal of the ab initio electronic structure calculations performed in this study was to obtain highly accurate *12*, *20*–*22*). In the case of the ground-state PEC, we combined the numerically exact description of the valence electron correlation effects provided by full CI with the high-level description of subvalence correlations involving all electrons but the 1s shells of Mg atoms obtained using CCSDT (*34*, *35*). Thus, the _{2} reported in this work was obtained by adopting the composite scheme

The first term on the right-hand side of Eq. 1 denotes the total electronic energy obtained in the full CCSDT calculations correlating all electrons other than the 1s shells of the Mg monomers and using the aug-cc-pwCVQZ basis set developed in (*40*), abbreviated in this section and in the Supplementary Materials as AwCQZ. The second and third terms on the right-hand side of Eq. 1, which represent the difference between the frozen-core full CI and CCSDT energies obtained using the aug-cc-pV(Q + *d*)Z basis of (*40*), abbreviated in this section and in the Supplementary Materials as A(Q + *d*)Z, correct the nearly all-electron CCSDT/AwCQZ energy for the valence correlation effects beyond CCSDT. The A(Q + *d*)Z and AwCQZ basis sets were taken from the Peterson group’s website (*41*). We used these bases rather than their standard aug-cc-pV*n*Z and aug-cc-pCV*n*Z counterparts, since it has been demonstrated that the aug-cc-pV(*n* + *d*)Z and aug-cc-pwCV*n*Z basis set families, including A(Q + *d*)Z and AwCQZ, accelerate the convergence of bond lengths, dissociation energies, and spectroscopic properties of magnesium compounds (*26*, *40*). The aug-cc-pV(T + *d*)Z, aug-cc-pwCVTZ, and aug-cc-pwCV5Z bases (*40*), abbreviated in this section and in the Supplementary Materials as A(T + *d*)Z, AwCTZ, and AwC5Z, respectively, and used in the auxiliary calculations discussed in section S1 to demonstrate the convergence of our computational protocol with respect to the basis set size (see tables S1 and S2), were taken from the Peterson group’s website (*41*) as well.

As shown in section S1, the AwCQZ and A(Q + *d*)Z bases are large and rich enough to provide spectroscopic properties of the magnesium dimer that can be regarded as reasonably well converged with respect to the basis set size, to within ~0.1 to 2 cm^{−1} for the experimentally observed *v*″ ≤ 13 levels and ~3 to 5 cm^{−1} for the remaining high-lying vibrational states and *D*_{e} (see, e.g., table S2). Ideally, one would like to improve these results further by extrapolating, for example, the nearly all-electron CCSDT energetics in Eq. 1, which are responsible for the bulk of the many-electron correlation effects in Mg_{2}, to the complete basis set (CBS) limit. Unfortunately, a widely used two-point CBS extrapolation (*42*) based on the subvalence CCSDT/AwCTZ and CCSDT/AwCQZ data, which are the only CCSDT data of this type available to us, to determine the CBS counterpart of the first term on the right-hand side of Eq. 1 would not be reliable enough. As demonstrated in (*26*) and as elaborated on in section S1 (see table S2), a CBS extrapolation using the AwCTZ and AwCQZ basis sets worsens, instead of improving, the *D*_{e}, *r*_{e}, and vibrational term values of the magnesium dimer compared to the unextrapolated results using the AwCQZ basis. As shown in table S2, the CBS extrapolation using the AwCQZ and AwC5Z basis sets would be accurate enough, but the CCSDT/AwC5Z calculations for the magnesium dimer correlating all electrons but the 1s shells of Mg atoms are prohibitively expensive. One could try to address this concern by replacing CCSDT in Eq. 1 by the more affordable CCSD(T) approach (*32*), resulting in*d*)Z basis sets rather than the poor-quality CBS extrapolation from the CCSDT/AwCTZ and CCSDT/AwCQZ information.

In principle, one could extend the above composite scheme, given by Eq. 1, to the electronically excited *37*, *38*), but the nearly all-electron full EOMCCSDT calculations using the large AwCQZ basis set turned out to be prohibitively expensive for us. To address this problem, we resorted to one of the CR-EOMCCSD(T) approximations to EOMCCSDT, namely, CR-EOMCCSD(T),IA (*39*), which is capable of providing highly accurate excited-state PECs of near-EOMCCSDT quality at the small fraction of the cost. Thus, our composite scheme for the calculations of the *d*)Z basis. Before deciding on the use of CR-EOMCCSD(T),IA, we tested other CR-EOMCC schemes (*43*) by comparing the resulting *G*(*v*′, *J*′) values with the available experimentally derived data reported in (*21*, *44*). Although all of these schemes worked well, the computational protocol defined by Eq. 3, with the CR-EOMCCSD(T),IA approach serving as a baseline method, turned out to produce the smallest maximum unsigned errors and RMSD values relative to experiment.

While the *20*, *21*), one might wonder whether the neglect of the post–Born-Oppenheimer and relativistic effects in our ab initio calculations could substantially affect our main conclusions. According to (*20*, *21*), the non-adiabatic Born-Oppenheimer corrections (BOCs) for the *21*), may have to be accounted for, but, based on the numerical data reported in (*21*), its magnitude is well within the uncertainty of the ab initio calculations reported in this work. According to (*30*), special relativity reduces the dissociation energy *D*_{e} characterizing the ^{−1}, i.e., the relativistic effects change the *D*_{e} by about 1%. However, our preliminary analysis using the modified version of the ab initio protocol adopted in the present work, in which the valence full CI and CCSDT calculations using the A(Q + *d*)Z basis set and the nearly all-electron CCSDT/AwCQZ computations are replaced by their scalar-relativistic counterparts employing the third-order Douglas-Kroll (DK) Hamiltonian (*45*, *46*) and the triple-ζ aug-cc-pV(T + *d*)Z-DK and aug-cc-pwCVTZ-DK bases (*40*), demonstrates that the number of bound vibrational states supported by the relativity-corrected *d*)Z and AwCTZ bases [the small negative differences between the relativity-corrected and nonrelativistic rotationless *G*(*v*″) values vary from <1 cm^{−1} or 0.8% for *v*″ = 0 to 2 to ~1% for the highest vibrational states near the corresponding dissociation thresholds]. Similar applies to the *D*_{e} value, estimated using the triple-ζ DK analog of the quadruple-ζ nonrelativistic computational protocol adopted in this work, is 0.2%, but the total number of bound vibrational states supported by the nonrelativistic and relativity-corrected potentials remains the same. The ab initio vibrational spectra corresponding to the *G*(*v*″ + 1) − *G*(*v*″) and *G*(*v*′ + 1) − *G*(*v*′) energy spacings do not exceed 0.4 cm^{−1} in the former case and 0.3 cm^{−1} in the case of the latter energy differences. Thus, while our preliminary findings regarding the small, but nonnegligible, effects of relativity need a thorough reexamination using both the larger basis sets, such as aug-cc-pV(Q + *d*)Z-DK and aug-cc-pwCVQZ-DK, and the various truncations in the DK Hamiltonian expansions, which may influence the calculated spectra too (*46*), and we will return to these issues in the future work, the *20*, *21*) and to comment on the corresponding rovibrational manifolds, especially for the ground electronic state.

All electronic structure calculations for Mg_{2} performed in this study, summarized in tables S3 to S5, were based on the tightly converged restricted Hartree-Fock (RHF) reference functions (the convergence criterion for the RHF density matrix was set up at 10^{−9}). The valence full CI calculations for the *47*), whereas the valence and subvalence CCSDT computations for the *48*). The valence and subvalence CR-EOMCCSD(T),IA calculations for the *39*), which take advantage of the underlying ground-state CC codes described in (*49*) and which are part of GAMESS as well. The GAMESS RHF-based CC routines (*49*) were also used to perform the CCSD(T) calculations needed to explore the basis set convergence and the viability (or the lack thereof) of the alternative to the CCSDT-based composite scheme given by Eq. 1, defined by Eq. 2 (see sections S1 and S2, especially table S2 and fig. S1). The convergence thresholds used in the post-RHF steps of the CC and EOMCC computations reported in this work were set up at 10^{−7} for the relevant excitation amplitudes and 10^{−7} hartree (0.02 cm^{−1}) for the corresponding electronic energies. The default GAMESS input options that were used to define our full CI calculations guaranteed energy convergence to 10^{−10} hartree.

The grid of Mg-Mg separations *r*, at which the electronic energies of the *r* values to determine the electronic transition dipole moment function *d*)Z basis set of (*40*) (see fig. S2 and table S5).

It is worth pointing out that our ab initio data points representing the *r* values are consistent with the expected long-range physics. One can see this, for example, by comparing our ab initio electronic energies for the *20*). Indeed, if we align the X-representation and our ab initio potentials such that the energies at *r* = 100.0 Å are identical (without this alignment, the X-representation and our ab initio energies at *r* = 100.0 Å calculated relative to the corresponding potential minima differ by 0.9 cm^{−1}), the differences between the two PECs in the *r* > 8.5 Å region, where the X-representation potential has the form ^{−1}, rapidly approaching zero as *r* increases (*r* = 100.0 Å is large enough to define the asymptotic region; for example, the difference between the X-representation energies at *r* = 30.0 and 100.0 Å is only 0.004 cm^{−1}; our ab initio calculations at the same two *r* values produce the numerically identical energy difference). In the *r* ≥ 20 Å region, where the X-representation energies are flat to within about 0.05 cm^{−1}, the X-representation and shifted ab initio energies, as described above, differ by less than 0.01 cm^{−1}.

### Calculations of rovibrational term values and rovibronic transitions

The rovibrational term values, including bound and quasi-bound states supported by our ab initio *50*) available in the LEVEL16 code (*51*) [LEVEL16 uses the Airy function approach described in (*52*) to locate quasi-bound states]. The widths and the tunneling lifetimes for predissociation by rotation characterizing the quasi-bound rovibrational states supported by the *51*) and implemented in LEVEL16, which requires numerical integrations between turning points in the classically allowed and classically forbidden regions of the relevant effective potentials including centrifugal barriers shown, for example, in fig. S3 [see (*51*) for further details]. To produce electronic energies *V*(*r*) on a dense grid of internuclear distances *r* with the step size of 0.001 Å, needed to perform the required numerical integrations and determine the corresponding equilibrium bond lengths *r*_{e} and dissociation energies *D*_{e}, we proceeded as follows. To obtain *V*(*r*) values every 0.001 Å in the *r* = 2.3 to 30.0 Å region, which excludes the innermost and outermost PEC parts defined by the 2.2 to 2.3 Å and 30.0 to 100.0 Å intervals, we used cubic splines available in LEVEL16, interpolating between pairs of nearest-neighbor *r* values used in the ab initio electronic structure calculations, starting from (2.3 Å, 2.4 Å) and ending up with (25.0 Å, 30.0 Å). To generate the equally densely spaced electronic energies in the innermost and outermost segments of each of the two PECs considered in this work, we resorted to analytical potential fits provided by the LEVEL16 code. In the case of the innermost parts of the *V*(*r*) = *A* + *Be*^{−Cr}, where parameters *A*, *B*, and *C* were determined by fitting the respective electronic energies calculated at *r* = 2.2, 2.3, and 2.4 Å. For the outermost, 30.0 to 100.0 Å, PEC segments, we adopted the appropriate long-range forms of the *C*_{6}, *C*_{8}, and *C*_{10} coefficients entering the former formula were obtained by fitting the *r* = 25.0, 30.0, and 100.0 Å to Eq. 4, in which *D*_{e} was defined as the relevant energy difference between *r* = 100.0 Å and *r*_{e}, with *r*_{e} representing the previously determined equilibrium internuclear separation in the ground electronic state. The six coefficients *C*_{3} through *C*_{8} entering the latter expression were obtained by fitting the *r* = 13.0, 15.0, 20.0, 25.0, 30.0, and 100.0 Å to Eq. 5, in which *D*_{e} was set as the energy difference between *r* = 100.0 Å and the corresponding *r*_{e}.

The quality of the potential fits generated by LEVEL16 is very high. We illustrate it here by summarizing the results of two of the several numerical tests that we carried out for the ground-state PEC. In one of the tests, we computed the electronic energy of the *r* = 3.893 Å, which is the equilibrium bond length determined by the potential fit *V*(*r*) produced by LEVEL16, using our ab initio quantum chemistry protocol defined by Eq. 1. The resulting energy, determined relative to the asymptotic value of the *r* = 100.0 Å, matched the value of *V*(*r*) at *r* = 3.893 Å obtained with LEVEL16 to within 0.0001 cm^{−1}. In another test, aimed at examining the ability of the interpolation scheme used by LEVEL16 to reproduce the ab initio energetics obtained using Eq. 1, we removed the electronic energies calculated at *r* = 4.6, 5.2, 5.8, and 6.4 Å, which is the region of the *V*(*r*) changes its curvature, and regenerated the potential fit using the remaining ab initio points. The new potential fit, based on fewer ab initio energies, reproduced the equilibrium bond length *r*_{e} = 3.893 Å resulting from the original fit, constructed using all values of *r* in our grid, to within 0.001 Å (which is the step for numerical integration in LEVEL16). The mean unsigned error characterizing the *r* = 4.6, 5.2, 5.8, and 6.4 Å resulting from the new fit relative to their ab initio values obtained using Eq. 1 was 0.05 cm^{−1}. On the basis of these and similar analyses, including the

We also used LEVEL16 to determine the rovibrational term values characterizing the experimentally derived analytical X-representation potential developed in (*20*), which we used to assess the accuracy of our ab initio–determined *r* adopted in our ab initio work. We then followed the same numerical procedure as described above for the

Last, but not least, we used LEVEL16 to compute the line positions of all allowed *v*′, *J*′) → *v*″, *J*″) rovibronic transitions and, with the help of our ab initio transition dipole moment function *v*′, *J*′) → *v*″, *J*″) transitions with the LIF data reported in (*20*, *21*) was a uniform downward shift of the entire ^{−1}, needed to match the experimentally determined adiabatic electronic gap *T*_{e} of 26,068.9 cm^{−1} (*21*). Other than that, all of the calculated spectroscopic properties, including the *D*_{e}, *r*_{e}, and rovibrational term values corresponding to the *v*′, *J*′) → *v*″, *J*″) transitions reported in this study, rely on the raw ab initio data compiled in tables S3 to S5 and data files S1 and S2, combined with the LEVEL16 processing, as described above. To produce Fig. 3, we superimposed our theoretical *v*′ = 3, *J*′ = 11) → *v*″, *J*″ = 10,12) LIF spectrum on top of the experimental one reported in figure 3 of (*20*). The theoretical line intensities shown in Fig. 3 were normalized such that the tallest peaks in the calculated and experimental LIF spectra corresponding to the *v*″ = 5 P12 line representing the *v*′ = 3, *J*′ = 11) → *v*″ = 5, *J*″ = 12) transition match.

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/6/14/eaay4058/DC1

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**

**Funding:**This work has been supported by the Chemical Sciences, Geosciences and Biosciences Division, Office of Basic Energy Sciences, Office of Science, U.S. Department of Energy (grant no. DE-FG02-01ER15228) and National Science Foundation (grant no. CHE-1763371).

**Author contributions:**I.M. and P.P. conceived the project. S.H.Y. and I.M. performed the electronic structure calculations under P.P.’s supervision. S.H.Y. carried out the rovibrational analysis and computed the rovibronic spectrum, which was subsequently used to simulate the relevant LIF spectra. All authors contributed equally and extensively to the discussion and interpretation of the data. All authors worked closely together on preparing the manuscript. P.P. supervised the entire project.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**All data needed to evaluate the conclusions of this study are included in the paper and/or the Supplementary Materials. Any additional information related to this paper may be requested from the authors.

- Copyright © 2020 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).