Calculation of Some Low-Lying Electronic Excitations of Barium Monofluoride Using the Eom-CC3 Method with an Effective Core Potential Approach

I. Introduction

The search for the electron’s electric dipole moment, as well as laser cooling with the goal of optical trapping has triggered renewed interest in polar diatomic molecules, such as barium and strontium monofluoride (BaF, SrF). On the theoretical side this interest is reflected in large-scale computational efforts: a comprehensive study of the low-lying excitations was performed using a relativistic Fock-space coupled cluster (CC) approach [6]. These calculations led to a much improved agreement with experimental data, in particular for the excitation energies $\Delta T_e$, i.e., the differences in potential energy curves at the respective minima. Agreement was found at the level of $100c{m}^{-1}$ for SrF and BaF (and larger for CaF), which represented a dramatic improvement over previous theoretical results. Subsequently, a fully relativistic treatment, including quantum electrodynamic corrections and basis set extrapolation led to an order-of-magnitude improvement for BaF [7], and agreement with experimental excitation energies $T_e$ at the level of $10c{m}^{-1}$. The all-electron Fock-space CC method has been shown recently to account for a number of properties, such as ionization energies [8], and hyperfine constants [9].

Since all-electron calculations are computationally challenging, particularly in the relativistic frameworks, it is of interest to explore whether effective core potentials (ECP) and corresponding basis sets can be combined with high-order CC methods to reach the $100c{m}^{-1}$ level precision with reduced computational effort. Our interest in this regard is the investigation of these molecules within a cryogenic matrix environment, such as argon or neon [3,10,11]. To investigate BaF in a such an environment a more economical approach is required, and this represents the motivation for this study.

The layout of the paper is as follows. We begin in Sec. Section II with a short summary of the basic equations which lead to the spectroscopic constants. In Sec. Section III we demonstrate how well the EOM-CC3 method combined with an effective core potential and appropriate basis set works for the $X{ }^2{\Sigma }^{+}\to B{ }^2{\Sigma }^{+}$ transition in SrF, which is free of spin-orbit interactions. We then apply the EOM-CC3, EOM-CCSD and state specific $\Delta -CC$ methods to the low-lying excitations in BaF. We draw some conclusions in Sec. Section IV.

II. Methodology

There are a number of coupled-cluster (CC) methods to compute electronic excitations. We have recently compared some of them for magnesium monofluoride and found that the EOM-CC3 method (EOM stands for equation of motion) when combined with the the aug-cc-pVnZ basis set family can reach agreement with experimentally derived electronic excitation energies $T_e$ [12]. The goal of the present work is to show that this methodology can work equally well for heavier systems (such as SrF or BaF) when reducing the all-electron problem to a (10+9)-electron problem by means of an effective core potential ECP [4,13]. We use the correlation consistent basis sets developed by Hill and Peterson [5].

To achieve high accuracy for electronic excitation energies of molecules the following two methods will be employed: EOM-CC3 with two basis set families, cc-pVnZ-PP and the augmented versions aug-cc-pVnZ-PP. For the former basis set family we were able to also apply the state-specific $\Delta$-CC methods which were reviewed in Ref. [14,15]. These methods can be economical computationally: one generates an unrestricted Hartree-Fock (HF) Slater determinant after orbital rotation (switching highest occupied and some low unoccupied orbital) to obtain a $\Delta$-SCF solution; then one follows up with a CC calculation for this state. The problem with this method is that for some large basis sets the HF calculation collapses back to the ground state. The computations in the $\Delta$-CC step for the excited state usually require many more iterations than for the ground state. An interesting feature, however, is the option to treat triple excitations perturbatively. i.e., there is a $\Delta$-CCSD(T) method for which there is no EOM counterpart.

One question we are interested in is how big differences are between the state-specific vs EOM methods, while using the cc-pVnZ-PP basis at the $n=5$ (5Z) level. The calculations were performed with the CCEOM package within Psi4 [16]. It allows for orbital rotations in the first irreducible representation, and we were able to carry out the necessary rotations by using different symmetries: $C_{2v}$, $C_s$, and $C_1$. Our focus is on the BaF excitations from $X{ }^2{\Sigma }^{+}$ to $A^{\prime }{ }^2\Delta$, $A{ }^2\Pi$, and $B{ }^2{\Sigma }^{+}$ states.

To obtain the spectroscopic parameters one has to calculate potential energy curves over some range of molecular separations $R_i$. A simpler task is to choose a fixed separation $R=R_e$ (e.g. the $X{ }^2{\Sigma }^{+}$ value of $R_e$ known from experiment), and to compute vertical excitation energies $\Delta T$. The $\Delta T$ values overestimate the excitation energies $\Delta T_e$ which are obtained from the potential energy curves for both ground and excited states using the respective state-specific equilibrium separations $R_e$.

The potential energy curves lead to the spectroscopic parameters ${\omega }_e$ and ${\omega }_e{x}_e$, which follow from the nuclear vibrational excitation energies in accord with [17]

$$E_v=T_e+{\omega }_e(v+\frac{1}{2})-{\omega }_e{x}_e{(v+\frac{1}{2})}^2.$$

(1)

A practical procedure we follow is to compute data for the potential energy curves for a set of values $R_i$, which span the region of R-values such that one obtains a realistic energy spectrum for low-lying nuclear vibrations of the molecule, e.g., $v=0,1,2,3$. A spline fit through the data allows to determine $R_e$, and the value of the minimum potential energy $T_e$.

It is convenient to fit these data to a Morse potential,

$$V\left(R\right)=T_e+D{(1-exp[-\beta (R-R_e)])}^2.$$

(2)

The two additional constants $D,\beta$ can then be found from the fit. Here D is just a parameter (not necessarily close to the dissociation energy), as the purpose of fitting to a Morse potential is to obtain a wider range in R to accommodate the tails of the vibrational states when solving the nuclear Schrödinger equation. This procedure avoids the computation at a larger grid of values $R_i$. Alternatively, one can carry out direct fits of the original data to the Morse potential eq.(2) to find the four parameters.

Fitting the vibrational spectrum for $v=0,1,2,3$ to eq. (1) allows one to find the parameters ${\omega }_e$ and ${\omega }_e{x}_e$, or a more complete set including ${\omega }_e{y}_e$ in case one uses not a Morse potential but a more general fit of the calculated potential energy points at the given set $\left\{R_i\right\}$. For the Morse potential the values of ${\omega }_e$ and ${\omega }_e{x}_e$ are related directly to the potential parameters.

One can also deduce the parameters from the experimentally determined values of the spectroscopic constants. In atomic units one has

$$D={\displaystyle \frac{{\omega }_e^2}{4{\omega }_e{x}_e}}\beta =\sqrt{{\displaystyle \frac{M}{2D}}}{\omega }_e,$$

(3)

where M is the reduced mass of the diatomic molecule.

In terms of expected accuracy for the excitation energies the goal is to be within $100c{m}^{-1}$ of experimentally deduced values. Is this percent-level (or better) goal achievable? The EOM-CC3 method is claimed to be good to $40meV$ level ($\approx 300c{m}^{-1}$) [18], in general. The reason why we can expect to do better than that for our cases is that the excitations are predominantly of single-electron type. In the case of MgF it was shown that the method when coupled with the augmented correlation consistent basis sets can deliver per mille accuracy for transition energies [12].

The results were obtained on multi-core workstations (e.g., Apple MacBook Pro with 8 cores, but during most of the computations the CPU was loaded only at the 40% level). Excited-state calculation using the EOM-CC3 method for a single root took of the order of 10-12 hours with the aug-cc-pV5Z-PP basis and about half the time for cc-pV5Z-PP.

III. Results

A. SrF($X{ }^2{\Sigma }^{+}\to B{ }^2{\Sigma }^{+}$)

As a first case we look into an excitation of SrF for which the equilibrium distances of ground and excited state are very close to each other, namely $R=2.075$ and $R=2.08$ Å respectively. Vertical excitation energies are shown in Figure 1 as a function of R over a range of values that are needed to determine basic spectroscopic parameters. The variation of the excitation energies with R is large for the aug-cc-pVnZ basis with $n=3$, i.e., TZ, but less for the QZ ($n=4$) and 5Z ($n=5$) calculations.

The obtained values of $R_e,\Delta T_e,{\omega }_e,{\omega }_e{x}_e$ for both states are compared in Table 1 with values derived from experiments, and with the Fock-space (FS) CC calculation of Ref. [6]. The variation of $\Delta T_e$ with the basis set quality parameter $n=3,4,5$ is not systematic, and the validity of extrapolation to the complete basis set limit may not be working well, and is therefore omitted.

The spectroscopic parameters ${\omega }_e$ and ${\omega }_e{x}_e$ were determined in the following way: the potential energy data for both states were obtained on a grid of 16 equidistant R values spanning $1.93<=R<=2.26$; these data were fitted to a local spline using the Interpolation function in Mathematica; the minimum position $R_e$ and energy $T_e$ was determined using FindMinimum; after subtraction of $T_e$ the data were then fitted to either a Morse potential, or to a fourth-order polynomial including orders 2, 3 ,4 and centered on the determined value of $R_e$. Using the Morse and polynomial fits the Schrödinger equation for the nuclear motion was solved with NDEigensystem in Mathematica using the reduced mass for ${ }^{88}Sr{ }^{19}F$, which is the most abundant isotope, and which was selected in experiments [19,20,21].

For the 5Z basis the present results are seen to work very well with the excitation energy $\Delta T_e$ obtained within $20c{m}^{-1}$, which is an excellent result. The focus on the excitation to the SrF($B{ }^2{\Sigma }^{+}$) state in this section is motivated by the case of BaF, for which one can argue that effective and true excitation energies may differ by over $100c{m}^{-1}$ [23]. The problem of state mixing of one excited state with a vibrationally excited neighboring state becomes a problem for BaF due to the closeness of the $A^{\prime }{ }^2\Delta$ and $A{ }^2\Pi$ potential energy curves. The potential energy curves for the sequence of CaF, SrF and BaF molecules are shown in Figure 1 of Ref. [6]. In contrast to BaF for SrF there should not be a problem with comparing effective experimental and theoretical parameters for the $B{ }^2{\Sigma }^{+}$ state.

A comment on the comparison with the calculations of Ref. [6] is perhaps in order. The scope of that work, of course, is substantially wider: the relativistic framework with an effective hamiltonian for the molecular mean-field approach is reduced to a two-component formalism (hence it is designated as the X2C-FSCC method). It includes part of the Breit interaction, and a basis set was constructed in a careful way with the doubly augmented d-aug-pVQZ relativistic basis of Dyall [25] expanded manually. This effort was required due to the treatment of the all-electron problem, and resulted in accurate spin-orbit splittings (which is not relevant for the $B{ }^2{\Sigma }^{+}$ state). In contrast, the present work describes relativistic effects by a pseudopotential for the inner 28 electrons of Sr (ECP28 based on multi-configuration Dirac-Fock and Breit interaction, cf. Ref. [13]), and uses a correlation-consistent basis set sequence designed for the atom [5]. Within the EOM-CC3 method this much simpler approach yields an improved excitation energy, and somewhat better values for the equilibrium distances.

B. BaF($X{ }^2{\Sigma }^{+}\to A^{\prime }{ }^2\Delta$, $X{ }^2{\Sigma }^{+}\to A{ }^2\Pi$, $X{ }^2{\Sigma }^{+}\to B{ }^2{\Sigma }^{+}$)

In order to compare ECP based calculations which are non-relativistic in nature except for taking into account the contributions from core electrons of the heavy atom via the pseudopotential we need to look at spectroscopic results where the spin-orbit splittings are removed. A detailed analysis of BaF low-lying electronic excitations led to the discovery of the lowest excitation, namely the $A^{\prime }{ }^2\Delta$ doublet [26]. The paper gives the excitation energies $\Delta T_e$ and spectroscopic parameters for the levels considered here, but experiment naturally includes the spin-orbit effect associated with $J=\frac{3}{2},\frac{5}{2}$ for $A^{\prime }{ }^2\Delta$, and $J=\frac{1}{2},\frac{3}{2}$ for $A{ }^2\Pi$. A simple approach to obtain reference values for comparison with the present results is to average the spin-orbit splittings. These values are shown in Table 2.

In a follow-up paper [27] a detailed analysis was reported for the vibration-rotation bands $v=0,1,2$, and effective constants were determined. These were derived from a model and resulted in the determination of spin-orbit coupling constants. The study was then complemented by a deperturbation approach [23], based on the notion that the excited states are predominantly 5d states. Another set of effective excitation energies was determined. The deperturbation concerns the fact that the $v=2$ vibrational excitation of the $A^{\prime }{ }^2\Delta$ state becomes close to the $v=0$$A{ }^2\Pi$ state. Due to spin-orbit couplings the $B{ }^2{\Sigma }^{+}$ state is strongly affected, and the excitation energies $\Delta T_e$ derived from experimental data may not be exactly the quantities expected from the theoretical potential energy curves.

More complete measurements analyzing higher vibrational states are reported in Ref. [28], and the spectroscopic parameters are given in Table 2 as reference values. Recently measurements and analysis were performed not only for the most abundant isotope ${ }^{138}Ba{ }^{19}F$, but also for ${ }^{136}Ba{ }^{19}F$ [29].

The deperturbed values of $\Delta T_e$ given in Table 3 of Ref. [23] agree with the simple averaging procedure shown in Table 2 for the $A^{\prime }{ }^2\Delta$ state, are almost $20c{m}^{-1}$ higher for the $A{ }^2\Pi$ state, and about $120c{m}^{-1}$ lower for the $B{ }^2{\Sigma }^{+}$ state ($T_{\Sigma }=13944.5$). This should be kept in mind when looking at the comparison with theory.

On the theoretical side we note that the relativistic FSCC method of Ref. [6] gives higher excitation energies for the properly coupled states with some overestimations at the $150c{m}^{-1}$ level for $A^{\prime }{ }^2\Delta$ and $A{ }^2\Pi$. The spin-orbit splittings are obtained to reasonable accuracy. This cannot be said for the CASSCF+MRCI method with spin-orbit coupling calculated with perturbation theory [30], and our own experience while using Molpro [31] in an all-electron approach was similar in this respect.

How useful are calculated vertical excitation energies $\Delta T$ as approximations for $\Delta T_e$? For the $A{ }^2\Pi$ and $A^{\prime }{ }^2\Delta$ states this turns out to work, but less so for the $B{ }^2{\Sigma }^{+}$ state due to the variation in $\Delta T\left(R\right)$ as will be shown below in Figure 3. Our EOM-CC3 results for $R=2.16$ Å, with augmented basis for both Ba and F are shown in Table 2. For the excitation to the $B{ }^2{\Sigma }^{+}$ state the difference between $\Delta T_e$ and $\Delta T$ is appreciable.

Overall, the results are quite good, however with substantial room for improvement for excitation to the $A^{\prime }{ }^2\Delta$ state due to the overestimation by over $200c{m}^{-1}$. For the $B{ }^2{\Sigma }^{+}$ state the result for $\Delta T_e$ supports the notion that theory should be compared to the value from the deperturbation analysis [23].

In order to assess the quality of the EOM-CC3 results we show in Figure 2 the calculated data with $T_e$ subtracted for the $X{ }^2{\Sigma }^{+}$ and $B{ }^2{\Sigma }^{+}$ states. Also shown are fits to the data using the Morse potential eq. (2). The data points for $R=2.01$ Å were removed from the fit, and disagree with the Morse shape for both states. The reason for this deviation may be associated with basis set problems (superposition error) for small values of R. The dashed magenta curves are the potentials derived from the experimental parameters $R_e,{\omega }_e,{\omega }_e{x}_e$ shown in Table 2 using eq. (3).

The agreement is not perfect: results from R-values to the left of the minima may be affected by the mentioned problem for the first data point, particularly in the case of the $X{ }^2{\Sigma }^{+}$ state. The values for ${\omega }_e{x}_e$ depend somewhat on the choice of data points used for the fits, and, thus, their values shown in Table 4 are, at best, estimates. The values for $R_e$ should be accurate as stated, and those for ${\omega }_e$ are deemed to be accurate to 2.5 digits.

In Table 5 we present results from different CC methods and two basis set combinations for the vertical excitation energies at $R=2.16$ Å. We focus on the comparison of basis sets aug-cc-pV5Z-PP and cc-pV5Z-PP for the barium atom combined with aug-cc-pV5Z for fluorine, since the latter combination allowed us to apply state-specific CC methods in addition to the EOM approach.

We can summarize the data as follows: the best calculations are given in the second row for EOM-CC3 with augmentation on both atoms. Removing augmentation on Ba results in literally no change for the excitation to $A{ }^2\Pi$, a small increase for the $A^{\prime }{ }^2\Delta$ state, and a $45c{m}^{-1}$ increase for the $B{ }^2{\Sigma }^{+}$ state. Interestingly, the excitation energies from the $\Delta$-CC3 method are higher than the EOM-CC3 results with the same basis on the order of $40-50c{m}^{-1}$. The EOM-CCSD calculations yield systematically higher excitation energies, and are, thus, deemed less useful. The $\Delta$-CCSD results, on the other hand are lower for the $B{ }^2{\Sigma }^{+}$ state. The $\Delta$-CCSD(T) method gives results close to those from the $\Delta$-CC3 method. In Table A1 in the Appendix the energies are given for completeness. One can observe that the actual energies for different methods disagree much more than the excitation energies.

Concerning the state-specific $\Delta$-CC results for excited states, which are obtained from higher roots than the ground state, we make some following observations. The $\Delta$-SCF solutions, which are required as an orbital basis to perform the CC steps are easy to find in $C_{2v}$ or $C_1$ symmetry for the $A{ }^2\Pi$ and $B{ }^2{\Sigma }^{+}$ states, by replacing the highest occupied with nearby lowest unoccupied natural orbitals (HONO vs $LUNO+j$) with $j=0,1$ using Psi4 terminology. For the $A^{\prime }{ }^2\Delta$ state we used $C_s$ symmetry and used the $LUNO+4$ to replace the HONO to obtain the $\Delta$-SCF Slater determinant. The CCSD(T) calculation generates the CCSD energy as an intermediate result, while the $\Delta$-CC3 calculations are done separately. The convergence properties of the CC correlation energy calculations were similar to ground-state calculations (on the order of 20 iterations) when working in $C_s$ or $C_{2v}$ symmetry.

In Figure 3 the excitation energies from EOM-CC3 with both basis sets are shown as a function of internuclear separation. The crosses show the data with augmentation on both atoms and are calculated on a fine grid of R-values (cf. Table A1). The curves show the strong dependence of $\Delta T$ vs R for the excitation to $B{ }^2{\Sigma }^{+}$, which is responsible for the difference between $\Delta T$ and $\Delta T_e$ for this state (cf. Table 2). The basis set combination which lacks augmentation on Ba gives nevertheless results of good quality, and may be preferred for more complicated situations (such as BaF within a cryogenic Ne matrix).

The $X{ }^2{\Sigma }^{+}$ ground-state energies for $R=2.16$ Å which are shown in the Appendix in Table A1 depend strongly on which CC method is used, and also show variation with the level of basis (aug/aug vs cc/aug for Ba/F). Nevertheless, the vertical excitation energies from the $\Delta$-CC methods differ from EOM-CC3 on the $100c{m}^{-1}$ scale for $\Delta$-CCSD and $30c{m}^{-1}$ for$\Delta$-CCSD(T). Basis augmentation at Ba changes this excitation energy by $50c{m}^{-1}$.

At the level of full EOM-CC3 we find that only the $A^{\prime }{ }^2\Delta$ state excitation energy changes noticeably (it comes out higher by about $30c{m}^{-1}$, which is on top of an $\approx 200c{m}^{-1}$ overestimate). For the sake of computational economy one may recommend this approach. Neglecting the triples (CCSD over CCSD(T) or CC3) leads to $\approx 100c{m}^{-1}$ differences in the excitation energies.

Figure 3. EOM-CC3 vertical excitation energies $\Delta T$ in $c{m}^{-1}$ as a function of R (in Å) for the electronic transitions $X{ }^2{\Sigma }^{+}\to A^{\prime }{ }^2\Delta$ (green), $X{ }^2{\Sigma }^{+}\to A^{\prime }{ }^2\Pi$ (blue), and $X{ }^2{\Sigma }^{+}\to B{ }^2{\Sigma }^{+}$ (red). The curves represent spline fits to data obtained with the basis cc-pV5Z-PP for barium and aug-cc-pV5Z for fluorine. The crosses are data points obtained with both basis sets augmented, which lowers the values for the excitation energies by up to $50c{m}^{-1}$.

Figure 3. EOM-CC3 vertical excitation energies $\Delta T$ in $c{m}^{-1}$ as a function of R (in Å) for the electronic transitions $X{ }^2{\Sigma }^{+}\to A^{\prime }{ }^2\Delta$ (green), $X{ }^2{\Sigma }^{+}\to A^{\prime }{ }^2\Pi$ (blue), and $X{ }^2{\Sigma }^{+}\to B{ }^2{\Sigma }^{+}$ (red). The curves represent spline fits to data obtained with the basis cc-pV5Z-PP for barium and aug-cc-pV5Z for fluorine. The crosses are data points obtained with both basis sets augmented, which lowers the values for the excitation energies by up to $50c{m}^{-1}$.

IV. Conclusions

EOM-CC3 calculations for electronic excitations of SrF and BaF molecules are reported on the basis of treating these molecules as 19-electron systems, while representing the inner cores of Sr and Ba with pseudopotentials, which effectively removed 28 or 46 electrons for the two heavy atoms.

For SrF only the $X{ }^2{\Sigma }^{+}\to B{ }^2{\Sigma }^{+}$ transition was investigated. The difference between the equilibrium nuclear separations $R_e$ for the two states is small. It was shown how the vertical excitation energy changes with separation R as a function of basis set parameter n. For the QZ and 5Z basis (aug-cc-pVnZ-PP/aug-cc-pVnZ for Ba/F) this variation is much less pronounced than for TZ. The resulting value for $\Delta T_e$ at the 5Z basis level is $20c{m}^{-1}$ below the experimentally determined value, which is substantially better than other calculations. The variation of $\Delta T$ for fixed R as a function of basis set quality parameter n with $n=3,4,5$ is not systematic, and therefore the extrapolation to the complete basis set limit for $\Delta T_e$ appears to be less valid than for the ground-state energy.

For BaF calculations were carried out for excitations to three closely related states, namely $A^{\prime }{ }^2\Delta$, $A{ }^2\Pi$, and $B{ }^2{\Sigma }^{+}$ using basis sets at the 5Z level. The $\Delta T_e$ values obtained for the $A^{\prime }{ }^2\Delta$ state overestimate the experimentally derived value by over $200c{m}^{-1}$, while the results for $A{ }^2\Pi$ and for $B{ }^2{\Sigma }^{+}$ agree to within $30c{m}^{-1}$ of the spectroscopic value obtained from a deperturbation analysis [23].

The finding of the present study for the $B{ }^2{\Sigma }^{+}$ excitation energy, namely that much better agreement is found with the result from a deperturbation analysis [23] than with the effective parameters derived directly from the experimental data of $v=0,1,2$ vibrational bands [27] is hard to reconcile with the theoretical analysis of Ref. [7], where agreement was found with the effective spectroscopic constant $\Delta T_e$ as quoted in Ref. [26].