An Intruder-Free Fock Space Coupled-Cluster Study of the Potential Energy Curves of LiMg+ within the (2,0) Sector

The potential energy curves (PECs) and spectroscopic constants of the ground and excited states of a LiMg+ molecular cation were investigated. We obtained accurate results for the fifteen lowest-lying states of the LiMg+ cation using the Intermediate Hamiltonian Fock Space Multireference Coupled Cluster (IH-FS-CC) method applied to the (2,0) sector. Relativistic corrections were accounted for using the third-order Douglas–Kroll method. In each instance, smooth PECs were successfully computed across the entire range of interatomic distances from equilibrium to the dissociation limit. The results are in good accordance with previous studies of this molecular cation. Notably, this study marks the first application of IH-FS-CC in investigating a mixed alkali and alkaline earth molecular cation, proving its usability in determining accurate PECs of such diatomics and their spectroscopic constants.


Introduction
In recent decades, researchers have shown interest in alkali and alkaline earth molecular dimers given their applications in studies conducted at ultralow temperatures, such as the controlled preparation of many-body entangled states [1] or the measurement of proton-to-electron mass ratio [2].Also, their respective molecular cations recently attracted the interest of experimental scientists in the context of ultracold studies.Various studies have demonstrated the usefulness of these systems: the KCa + cation was applied in the quantum simulation of solid-state physics [3], NaCa + cations can be used in quantum information processing as quantum logic gates, and RbBa + can be used to study strong-coupling polaronic effects [4,5].Furthermore, previous findings show that ion-atom sympathetic cooling in a magneto-optical trap is possible for NaCa + cations [6].The LiCa + cation can be applied as a high-spatial-resolution probe in ultracold chemical reactions [7], and the LiMg + cation can be used to measure the electron dipole moment [8].Moreover, RbCa + and RbSr + cations can both be applied in studies of collisional quantum features such as scattering resonances or intermolecular effects [9,10].
A highly precise knowledge of potential energy curves (PECs) is essential for the evaluation of the properties of a system and for understanding collision and dissociation processes.Moreover, theoretically computed PECs provide necessary information for experimentalists, aiding them in adjusting their equipment to obtain ultracold species [11].However, calculating PECs using standard computational schemes can be challenging, even with relatively small systems like diatomic molecular cations.Difficulty arises when treating the homolytic dissociation of a single bond in which open-shell fragments are produced.Therefore, using the restricted Hartree-Fock scheme (RHF) for distances far from equilibrium is incorrect.Alternatively, unrestricted Hartree-Fock (UHF) or restricted openshell Hartree-Fock (ROHF) methods are required, which are known for their convergence problems with the HF and post-HF equations.Theoretical chemistry has proposed several solutions to address this issue.A frequently employed method in studies of excited states (EEs) is the Equation-of-Motion Coupled-Cluster (EE-EOM-CC) method [12][13][14][15][16][17].Unfortunately, this method is not size-extensive, which is the reason for its limited use in studies of PECs [18,19].A more commonly used approach is the multireference configuration interaction (MRCI) method, which is size-extensive only in its full CI (FCI) variant, making it impossible to carry out in a reasonable time period for systems with large numbers of electrons.As a workaround, researchers often opt to perform calculations exclusively for valence electrons while freezing the core, impacting the accuracy of the results.To overcome this limitation, some studies use pseudopotential or ECP (Effective Core Potential) approximations [20][21][22].In this approach, electron correlation is accounted for only for valence electrons, and it is usually reduced to the CISD (S-singles; D-doubles) method, which is to the FCI for the two-electron case.This approach also requires considering some additional parameters representing the potential of core electrons, and they cannot be used universally for any chosen system.
In this study, we focus on computing the accurate potential energy curves and spectroscopic constants of a molecular cation, LiMg + .We decided to utilize an alternative approach using the IH-FS-CCSD method (the Intermediate Hamiltonian Fock Space Multireference CC method with singles and doubles) [23], obtaining results in the (2,0) sector, which is explained further in the Methods section of this manuscript.This method allows one to choose a doubly ionized system as a reference, dissociating it into closed-shell fragments.This makes it possible to use the RHF as a reference function for any internuclear distance.Then, calculations using the double-electron attachment (DEA) formalism are performed to produce energies of the ground and excited states of the studied system.In the case of the system studied here in the form of the LiMg + molecular cation, the reference function was determined for the triple-positive ion LiMg 3+ : then, utilizing the DEA formalism: we obtain the energies of the ground and excited states of the LiMg + .This approach allows smooth PECs to be obtained for the whole range of internuclear distances.IH-FS-CCSD(2,0) is also a strictly size-extensive method-this property ensures that the energies of the electronic states of the system converge at an infinite distance to the sum of its atomic values.
The IH-FS-CCSD(2,0) method [23] was previously proven successful in obtaining PECs and the spectroscopic constants of diatomic molecules composed of alkali metals in several studies, as demonstrated in recent papers involving NaLi [24] or LiRb [25] molecules.However, this marks the inaugural application of this computational approach to study a molecular cation composed of alkali and alkaline earth metals, specifically the LiMg + molecular cation.Our objective was to define a suitable methodology for the calculation of such species and present benchmark accuracy results.
The next section of this work presents computational details and the results obtained, along with a discussion.The Methods section describes the theoretical background of the IH-FS-CCSD(2,0) method, and the final section offers conclusion, new insights, and future perspectives.

Background
The LiMg + system was previously studied with respect to its spectroscopic properties, both theoretically and experimentally, in only a few papers.The oldest theoretical works concerned only the spectroscopic constants of the ground state of the molecular cation [26][27][28].More recent papers published in the last decade employed various methods to study the spectroscopic constants of the ground and excited states of LiMg + , including the MRCI calculations of Gao and Gao [29], the pseudopotential study of ElOualhazi and Berriche [8], the calculations of Fedorov et al. using CC and MRCI methods [30], the EE-EOM-CC study of Bala et al. [18], and the CC calculations of Śmiałkowski and Tomza [11].In terms of experimental studies, there is a single paper by Persinger et al. in which the spectroscopic values of only the ground state of LiMg + are given since the paper was focused on the neutral LiMg molecule [31].

Computational Details
Calculations were carried out using the IH-FS-CCSD(2,0) method [23], which is explained in detail in the Methods section of this manuscript.The main results were computed using the uncontracted ANO-RCC [32] basis set with additional diffuse functions, which we called unANO-RCC+.For each atom, we have added six diffuse functions which are presented in Table 1.These exponents were determined through the even-tempered scheme [33]-the ratio between the subsequent exponents is equal to 0.35 for lithium and 0.4 for magnesium.This ensures the accurate ordering of atomic states.The final dimension of the unANO-RCC+ basis set is 242 spherical harmonic basis functions.The LiMg 3+ cation was chosen as the reference system for all double-electron attachment calculations, and all electrons were correlated.The reference function employed throughout this study was always the restricted Hartree-Fock function.To assess the impact of relativistic effects on the spectroscopic constants of the LiMg + cation, we adopted a two-step approach.Initially, the results for the LiMg + cation were computed using the uncontracted ANO-RCC basis set without the additional diffuse functions that are presented above-we called this basis unANO-RCC.Afterward, scalar relativistic corrections were applied using the third order Douglass-Kroll [34] method, and the calculations were repeated using the standard ANO-RCC basis set.The differences between obtained relativistic and non-relativistic values were computed and added to the spectroscopic constants derived from the unANO-RCC+ basis set.

Atomic Energies at the Dissociation Limit
Size extensivity is a crucial property in PEC calculations-the energies of the electronic states of a system must converge to the sum of their atomic values at an infinite distance.IH-FS-CCSD(2,0) is a strictly size-extensive method.It is presented in Table 2, which displays values computed using the three basis sets that were utilized in this work.The first column depicts a given dissociation limit, while the next columns depict energy values for Li or Li + and Mg or Mg + , respectively.The second-to-last column displays the sum of these energies, while the last column presents the energy at the corresponding dissociation limit obtained using the IH-FS-CCSD(2,0) method.The energies of the Li + cation were computed using the CCSD method (no electrons attached) for the Li atom; for the Mg + cation, we used IH-FS-CCSD(1,0) equivalent to EA-EOM-CCSD (one electron attached); and for the Mg atom, the IH-FS-CCSD(2,0) method was utilized (two electrons attached).The results are equal, proving the size extensivity of the method.We also compared the obtained excitation energies of the lithium and magnesium atoms with available experimental data [38], as shown in Table 3.In the columns, we depict the four lowest excitation states of Li and the same for Mg; in the rows, there are values for the three basis sets used in this work.The agreement between the experimental values is dependent on the basis set.For lithium, the ANO-RCC basis set performs accurately only for the first two excited states, altering the order of the third and fourth.The unANO-RCC+ and unANO-RCC basis sets are more accurate, with differences ranging from less than 0.001 eV to 0.003 eV and from 0.001 eV to 0.084 eV, respectively.For magnesium, ANO-RCC is more accurate, with discrepancies ranging from 0.011 eV to 0.042 eV.The differences for unANO-RCC+ and unANO-RCC are between 0.015 eV and 0.025 eV and between 0.002 eV and 0.028 eV, respectively.Importantly, only the first two terms of lithium are a part of the dissociation limit presented in this work; therefore, we were able to use the ANO-RCC basis set.In summary, our method closely reproduces the experimental excitation energies.

Potential Energy Curves
The potential energy curves of the fifteen lowest-lying electronic states of LiMg + were calculated using the IH-FS-CCSD(2,0) method and the unANO-RCC+ basis set, and they are presented in Figure 1.The curves are distinguished by eight different colors corresponding to eight dissociation limits and four different point types depending on the multiplicity and symmetry of their states.The interatomic distances presented in Figure 1 are limited to 25 Å, but the tabulated total energies of each state up to 200 Å are provided in the Supplementary Materials.
Among the fifteen PECs presented in this work, the 4 1 Σ + state is the only one characterized as having two minima.As expected, the shapes of higher-lying curves are distorted, with noticeable undulations, and they do not represent regular Morse-like potentials due to the phenomenon of avoided crossing between adjacent adiabatic potential energy curves of the same symmetry and multiplicity.Its presence can be explained by the occurrence of the interaction and charge transfer process between the electronic states of Li-Mg + and Li + -Mg.To better illustrate these interactions, we present specific PECs in Figure 2. The exact positions of the avoided crossing interactions are shown in Table 4, and they are compared with positions presented in [8].The agreement is very good, with differences of less than 2%.We also depicted PECs obtained using the IH-FS-CCSD(2,0) DK3/ANO-RCC method in Figure 3.The tabulated total energies of each state up to 300 Å are also provided in the Supplementary Material.These PECs are shown purely for comparison as they were used only while evaluating the influence of relativistic effects on the spectroscopic constants.Nonetheless, their shapes closely resemble those obtained using the unANO-RCC+ method (cf. Figure 1).Furthermore, the PECs presented in Figures 1 and 3 uniformly dissociate into the same limits, and their shapes align with those available in the literature [8].

Spectroscopic Constants
The following spectroscopic constants of the LiMg + molecular cation were obtained: equilibrium distances R e , well depths D e , adiabatic excitation energies T e , harmonic frequencies ω e , anharmonicity constants ω e x e , and equilibrium rotational constants B e .The results obtained using the IH-FS-CCSD(2,0)/unANO-RCC+ method are shown in Table 5.The relativistic contributions for each constant are given in parentheses.The magnitude of the impact of relativistic effects on the spectroscopic constants of the LiMg + cation varies for different states, but we observe that the contribution is smaller for lower-lying excited states than for higher-energy states.These relativistic contributions are smaller than for, e.g., LiRb [25], but they are still non-negligible, and they have some impact on the quality of the results.
Comparing our results for ω e and ω e x e with the paper by Persinger et al. [31], we can clearly see that our calculations reproduce the experimental values very well, with absolute differences of less than 1% and ∼7.7% for ω e and ω e x e , respectively.Moreover, our ω e value is the closest to the experimental one among all previous theoretical works.
Our results are also in good agreement with theoretical values obtained in previous works for most cases.The exceptions are values calculated for the 2 1 Σ + and 3 1 Σ + states in a paper by Bala et al. [18], in which they explicitly stated that their R e was largely different from that considered in [8], which makes it impossible to compare the values directly.
For the ground state, our R e value is the same as the CCSDT value obtained by Fedorov et al. [30].Also, our D e value is the closest to the result [11] calculated using the so-called "gold standard" CCSD(T) method [39].The IH-FS-CCSD(2,0) method is capable of providing very accurate results even if it does not consider the effect of triple excitations.a the method used in this work, which is IH-FS-CCSD(2,0)/unANO-RCC+ plus DK3 relativistic correction, as described in detail in the text.b the method used in [11] is CCSD(T)/aug-cc-pCVQZ. c the first method used in [30] is CCSDT/cc-pCVQZ.d the second method used in [30] is MRCI/cc-pCVQZ.e the method used in [18] is CCSD(T)/cc-pVQZ with a relativistic correction for the ground state and EE-EOM-CCSD/cc-pVQZ for excited states.f the methods used in [29] are the MRCI plus Davidson and third-order Douglas-Kroll-Hess relativistic correction using the aug-cc-pV5Z basis set.g the method used in [8] is based on pseudopotentials.The basis set for lithium was (9s8p5d/8s6p3d); for magnesium, it was (9s7p5d4f/7s7p4d4f).
The comparison with the paper by ElOualhazi and Berriche [8] is the most extensive as it is the only work comprising values for all of the excited states of LiMg + studied here.The closest R e values are observed for the 1 3 Σ + and 3 1 Σ + states, for which the absolute difference is 0.002 Å.The most significant error-less than 8%-is calculated for the 6 1 Σ + state.The good alignment of R e indicates a similar agreement in the derived B e values.
Also, we identified one repulsive state, i.e., the 2 3 Π state.The cited paper defines this state as having an extremely shallow D e of just 2 cm −1 ; thus, we believe it is probably a plateau.Moreover, we determined two states-4 1 Σ + and 2 1 Π-to have barriers, both of which retain the same dissociation limit.The heights of these barriers are 212 cm −1 and 91 cm −1 , respectively.These values correlate to the D e values reported in the cited paper.For the other states, our D e is closest in the case of the 1 1 Π state (an absolute difference of 20 cm −1 ) and farthest in the case of the 5 3 Σ + state.
The T e values show the best agreement throughout all states among the considered spectroscopic constants; the largest difference is less than 3% for the 5 3 Σ + state, and for most states, it is less than 1%.In nearly all cases, our ω e values are similar to the cited ones.The exceptions are the 1 1 Π state and the inner minimum of the 4 1 Σ + state, for which the absolute differences are larger than 10 cm −1 .As for ω e x e , the agreement is weak in most cases; this is due to the shapes of the PECs in large internuclear distances, which are described well by the IH-FS-CCSD(2,0) method but significantly less well by pseudopotential-based methods.
In summary, the accuracy achieved is reliable in comparison with Ref. [8], though it varies depending upon the specific state and the spectroscopic constant in consideration.

Methods
Within the single-reference (SR) coupled-cluster theory, the Schrödinger equation is solved, assuming the exponential form of the wave function |Ψ⟩: where H is a Hamilton operator, E is the total energy of the system, |Φ 0 ⟩ is a reference (usually Hartree-Fock) function, and T is a coupled-cluster excitation operator.Within the current approximation of singles and doubles, it takes the form where the operator T 1 = ∑ ia t a i a † i is responsible for single excitation, i.e., moving an electron from an occupied level i to an unoccupied level a and, analogously, the T 2 = ∑ ijab t ab ij a † b † ji is responsible for double excitation, i.e., moving two electrons from occupied levels i,j to unoccupied levels a,b.In situations in which the reference function |Φ 0 ⟩ is of a multiconfigurational character (which is the case for potential energy curve calculations) it is necessary to generalize the SR approach to the multireference (MR) case, i.e., to replace the Schrödinger equation, Equation (3), with its Bloch equivalent: The superiority of the multireference approach relies on the fact that the exact energy, E, is obtained by solving the above equation or, equivalently, by the diagonalization of the effective Hamiltonian operator H e f f : within a small subspace (M) of the full configurational space (M 0 ), which is called a model space and is defined by the projection operator P, note that |Ψ 0 ⟩ is a component of the exact wave function |Ψ⟩ residing within the model space, i.e., |Ψ 0 ⟩ = P|Ψ⟩.Of course, the construction of the H e f f is not trivial, and it requires employing the wave operator Ω, reproducing, by definition, the exact wave function |Ψ⟩ by operating on the model function The form of the wave operator Ω depends on the assumed variant of the multireference theory.The characteristic feature of the Fock space formulation [41,42] used in the current work is a specific form of the model space which, by definition, can contain a variable number of electrons.Within the scheme adopted here, the valence Fock space is composed of three sectors defined by the following projection operators, respectively: P = P(0,0) + P(1,0) + P(2,0) ( The 0-valence sector corresponds to the reference function |Φ 0 ⟩, and its projection operator is equal to P(0,0) = |Φ 0 ⟩⟨Φ 0 |.The projectors for the remaining two sectors, one-valence P(1,0) and two-valence P(2,0) , are written as where |Φ α ⟩ represents the configuration with one additional electron placed on the valence level α, while |Φ αβ ⟩ represents the configuration with two additional electrons placed on the valence levels α and β.
For the (1,0) sector, we obtain and analogously, for the (2,0) sector, we have Owing to the killer condition, Equation ( 14), we eliminate from the exponential in Equation ( 16) the operators S (1,0) and S (2,0) , and, similarly, from Equation ( 17), we eliminate the operator S (2,0) .Restricting our model to singles and doubles (cf.Equation ( 5)), we have the following for the one-valence sector: and for the two-valence one: After solving the FS equations for the (0,0) sector, we may further simplify the equations for the (1,0) and (2,0) sectors by replacing the H operator with a similarity transformed Hamiltonian H defined as H = e −S (0,0) He S (0,0) = e −T He T Equations ( 17) and ( 18) now take form Q(2,0) He S (1,0) +S (2,0) P(2,0) = Q(2,0) e S (1,0) +S (2,0) H e f f P(2,0) The FS-CC method described here is usually referred to as a standard or H e f f -based scheme.In this approach, at each stage, one can encounter real problems with convergence due to the appearance of intruder states [43,44], i.e., determinants belonging to the orthogonal space with energies very close to those of the model determinants.Hence, the choice of the active space is, in fact, limited to small spaces.The larger the active space, the more difficult the solution is.
The intermediate Hamiltonian formalism [45][46][47] provides a way to avoid the troublesome iterative procedure of the FS problems and provides eigenvalues identical to those from the standard effective Hamiltonian-based method.Thus, the IH scheme can replace the iterative solution of the FS equations with a diagonalization of the properly constructed matrix; thanks to this, all limitations originating from the standard formulation disappear.
The main idea of the IH approach relies on selecting a part of the orthogonal space M ⊥ as an intermediate space M I connected by the PI projector.In our case, the intermediate space is spanned by all determinants which are obtained by the direct action of S (2,0) 2 on M (2,0) .For the IH-FS-CCSD(2,0) case [23], we have the subspaces M (2,0) , M Q(2,0) = P(2,0) In this approach, we introduce the following wave operators: which, at the IH-FS-CCSD(2,0) level take the forms (2,0) 2 (1,0) 2 Thus, in summary, the Fock space coupled-cluster approach has two realizations in the (2,0) sector.The first one, which we denote FS-CCSD(2,0), is based on solving Equations ( 22) and ( 23), and it is impractical to use due to the intruder problem.The second one, which we denote IH-FS-CCSD(2,0), relies on the construction and diagonalization of the H (2,0) I matrix.Note that both approaches give identical solutions as far as the energies of the electronic states are concerned, and the method we use to compute the PECs and spectroscopic constants is the IH-FS-CCSD(2,0) method.
One should also mention that the Fock space approach at the one-valence level, i.e., the (1,0) and (0,1) sectors, is equivalent to the EOM-CC scheme applied to electron affinity and ionization potential cases, respectively.This equivalence means that the eigenvalues obtained via both approaches are identical, while the eigenvectors can be obtained from each other via a simple transformation.The EOM level, i.e., the final stage of the whole sequence of computational steps, is quite straightforward and relies on the diagonalization of the H operator within the configurational subspace selected in accordance with the problem under investigation.For the (1,0) sector, we have the following: and i refers to the number of occupied one-particle levels.Thus, the key point here is the intermediate Hamiltonian formalism, which offers an easy way to replace the iterative solution of the standard Bloch equation, Equation (23), with the direct diagonalization of the IH matrix constructed on the basis of the amplitudes from the (1,0) sector and Equations ( 31) and (32), obtained via EOM method, see Equation (34).The final outcome of this is the elimination of the intruder state problem in the FS approach.Moreover, the size-extensive Fock space approach has a built-in capability to provide, by selecting a proper FS sector, the correlated results for an altered (with respect to the Hartree-Fock (HF) reference) number of electrons.For example, the FS(m, 0) sector produces results pertaining to a system with m electrons added to the HF function.Assuming that we take the neutral molecule as an HF reference, these results would correspond to the m-tuply negative anion.We may take advantage of this capability of the FS approach to choose-at the HF level-such a reference with a convenient property to dissociate into closed-shell fragments and to generate a smooth curve for the ground and excited states in the whole region of interatomic distances, keeping in mind that by selecting a proper sector of the Fock space, we will recover the original structure we want to study.

Conclusions
In this study, we performed accurate quantum-chemical calculations of the potential energy curves of a LiMg + molecular cation using the multireference coupled-cluster method, the IH-FS-CCSD(2,0) method.We overcame the challenge of properly describing the dissociation of closed-shell species into open-shell parts by employing the doubleelectron attachment formalism.Our calculations for the ground and excited states covered, for the first time, the entire range of interatomic distances using the RHF reference function, correlating all electrons.This approach also involves the Intermediate Hamiltonian formulation, which allows for the direct diagonalization of a suitably constructed matrix, mitigating complications associated with intruder states instead of finding an iterative solution, as in the standard Bloch equation.
For the first time, the IH-FS-CCSD(2,0) formalism was used to study the PECs and spectroscopic constants of a molecular cation composed of alkali and alkaline earth metals.We obtained smooth and accurate PECs of the fifteen lowest-lying states of the LiMg + cation approaching eight dissociation limits.Through this approach, we also obtained accurate spectroscopic constants.Among all previous theoretical studies, our ω e value of the ground state exhibits the closest agreement with the experimental value.Furthermore, for the ground state, we obtained a level of accuracy comparable to those of CCSD(T) and CCSDT studies without the need to include the effect of triples.For the excited states, our results indicate that most of our computed values align well with those from prior theoretical papers, and observed disparities could likely be attributed to the constraints of the EE-EOM-CCSD and pseudopotential-based methods.
Based on the quality of results obtained in various prior studies involving diatomic molecules comprising alkali metals [23][24][25]48,49] using the IH-FS-CCSD(2,0) method, we anticipate that the accuracy of presented the PECs and spectroscopic constants of LiMg + can be regarded as benchmark-level.
The present study aims to lay the groundwork for promising perspectives for future spectroscopic investigations of molecular cations of alkali and alkaline earth metals.It especially seeks to open a path for further theoretical and experimental studies on LiMg + , mainly in the context of research on ultracold species where accurate PECs and spec-troscopic constants are much needed.A wide range of potential applications of these molecular cations, e.g., in the studies of ultracold chemical reactions, quantum information processing, or measurements of fundamental physical constants, will be of great importance in the future.

Table 1 .
Additional diffuse functions for the Li and Mg atoms in the unANO-RCC+ basis set for the s, p, and d shells.

Table 2 .
Energies of electronic states at dissociation limit of LiMg + molecular cation compared to atomic energies in different basis sets.

Table 3 .
Excitation energies of Li and Mg atoms.All values are in eV.

Table 5 .
Spectroscopic constants of LiMg + .D e , T e , ω e , ω e x e , and B e are given in cm −1 ; R e is given in Å. DK3 relativistic corrections are given in parentheses.