Hyperfine Structure Constants on the Relativistic Coupled Cluster Level with Associated Uncertainties.

Accurate predictions of hyperfine structure (HFS) constants are important in many areas of chemistry and physics, from the determination of nuclear electric and magnetic moments to benchmarking of new theoretical methods. We present a detailed investigation of the performance of the relativistic coupled cluster method for calculating HFS constants within the finite-field scheme. The two selected test systems are 133Cs and 137BaF. Special attention has been paid to construct a theoretical uncertainty estimate based on investigations on basis set, electron correlation and relativistic effects. The largest contribution to the uncertainty estimate comes from higher order correlation contributions. Our conservative uncertainty estimate for the calculated HFS constants is ∼5.5%, while the actual deviation of our results from experimental values is <1% in all cases.

Introduction

The hyperfine structure (HFS) constants parametrize the interaction between the electronic and the nuclear electromagnetic moments. The HFS consequently provides important information about the nuclear as well as the electronic structure of atoms and molecules and can serve as a fingerprint of, for example, transition metal complexes, probed by electron paramagnetic resonance (EPR) spectroscopy,[1] or of atoms, ions, and small molecules in the field of atomic and molecular physics, investigated by optical or microwave spectroscopy. With the ever relentless progress in the field of atomic and molecular precision experiments, there is a growing need for both experimental and theoretical determination of the HFS. Accurate calculations of the HFS parameters can serve a direct as well as an indirect purpose, as will be elaborated in the following.

One example of a direct application of accurate theoretical HFS parameters is nuclear studies, where the calculated electronic properties (magnetic induction and electric field gradient) are used to extract the nuclear magnetic dipole and electric quadrupole moments of the heaviest or unstable atomic nuclei from the measured magnetic-dipole, A, and electric-quadrupole, B, HFS constants, respectively.[2,3] Another example is in the search for even better atomic clocks where the structure of the hyperfine levels must be known to great accuracy in order to make reliable predictions to guide new experiments.[4]

The calculated values of the HFS constants can be also used as a means to benchmark the employed theoretical method against existing experimental or higher level theoretical data. In order for a theoretical method to yield accurate predictions of the HFS constants, the electron distribution in the vicinity of the atomic nucleus in question must be properly described; comparison to experiment can thus give an indication of the quality of the employed wave function. Such applications can be considered to serve an indirect purpose.

Using HFS constants as benchmarks is particularly valuable when one is interested in a property that is sensitive to the interaction between electrons and nuclei and that cannot be obtained experimentally. One such example is the interpretation of the atomic parity nonconserving (PNC) measurements in Cs atoms, where theoretically determined PNC matrix elements are needed in order to extract the weak charge, i.e., the strength of the neutral weak interaction, from the measured transition amplitudes.[5,6] These matrix elements are sensitive to relativistic effects, which become important when the electrons are close to the atomic nucleus. Therefore, the accuracy of the calculated HFS constants (compared to experiment) serves as a good indication of the reliability of the predictions for the PNC matrix elements. In order to unambiguously test agreement with the Standard Model prediction of the weak charge, the uncertainty of the theoretical predictions needed to be smaller than 1%; such accuracy eventually was reached by several groups using calculated HFS constants as benchmark values.[7−12] Such system-specific sensitivity or enhancement factors are generally needed in the search for physics beyond the Standard Model with atoms and molecules.[13−15]

When accurate predictions of the HFS constants for heavy atoms or for molecules containing heavy elements are needed, special attention must be paid to two aspects: relativistic effects and electron correlation. In addition, it is desirable to use a method that allows reliable uncertainty estimates. In this study we investigate a scheme that meets these three requirements.

In the rest of this paper we will consider the magnetic-dipole HFS constant, which we will refer to as simply the HFS constant. We begin with an overview of the currently popular methods used in the calculations of this property.

As we are interested in high accuracy treatment of correlation and relativistic effects, we will limit this overview to methods that treat relativity beyond scalar relativistic effects and correlation beyond density functional theory (DFT). For an overview of nonrelativistic as well as DFT based methods we refer to the chapter by H. Bolvin and J. Autschbach.[16] For atoms, methods such as the multiconfigurational Dirac–Fock (MCDF),[17] Dirac-Hartree–Fock augmented by the many body perturbation theory (DHF + MBPT),[18,19] configuration interaction with MBPT (CI + MBPT),[20,21] all order correlation potential,[22] coupled cluster singles doubles with partial triples (SDpT)[8,23] and Fock-space coupled cluster (FSCC)[24] were shown to provide reliable results. For molecules, the situation becomes more complicated due to the lack of spherical symmetry and a limited number of implementations exist. These include the multireference configuration interaction (MR-CISD) method,[25] the restricted active space CI (RAS-CI) approach,[26] and the coupled cluster singles and doubles (CCSD) method.[26,27]

In this work we investigate the performance of the relativistic coupled cluster (CC) method for calculating the HFS constants of atoms and molecules. Where applicable, this approach provides the highest level of theory, while still being feasible for computations on the heaviest elements. In addition, the systematic construction of the CC method allows for a reliable uncertainty estimation. In this work we combine the CC approach with the well-known finite field scheme (also known as the finite difference method) to extract the HFS constants. This provides us with a straightforward way to calculate molecular properties as numerical derivatives.[28] The finite field approach is particularly useful in the framework of CC theory, since the formulation of expectation values is cumbersome due to the complicated form of the wave function. That said, several implementations exist for calculation of CC expectation values; the recent relativistic examples are the extended CC method (ECC),[26] the Z-vector CC method,[29] and analytic gradients approach.[30] An advantage of using the finite field method is that no truncation of the CC expansion is necessary (which is the case for the ECC method for example) and that it allows inclusion of the perturbative triple excitations without additional complications. A drawback of the finite field method is the increased computational cost. Furthermore, one has to pay special attention to the numerical stability.

The combination of the relativistic CC method and the finite field approach has previously been applied to various properties, such as dipole polarizabilities,[31] electric field gradients,[32−34] contact densities for calculating Mössbauer isomer shifts[35] and P- and P,T-odd relativistic enhancement factors.[36−38] The combination of the CC method and the finite field approach for calculating HFS constants has previously been used in a nonrelativistic framework,[39−41] but, to the best of our knowledge, the extension to a relativistic framework and application to systems with heavy atoms has not been demonstrated before. Here, we investigate the performance of this method and the effect of various computational parameters (e.g., basis set quality, active space size, treatment of higher order relativistic effects, and others) on the obtained results. Furthermore, we employ a straightforward and reliable scheme for assigning uncertainties of the calculated HFS constants.

Inspired by the examples mentioned above, we have chosen to apply our investigations to the HFS constants of the Cs atom and the BaF molecule. Due to the atomic PNC experiments, the HFS constant of Cs has been studied extensively and on high levels of theory, which makes it an ideal system for benchmark calculations. The BaF molecule is currently used in various experiments searching for physics beyond the Standard Model,[42−44] where theoretically determined enhancement factors are crucial for the interpretation of the measurements and the calculated HFS constants can provide an important indication of the theoretical uncertainty.

Theory

The magnetic hyperfine interaction between the electronic spin and the nuclear spin of the Mth nucleus is parametrized by the 3 × 3 hyperfine coupling tensor, A. It is usually defined through the effective spin Hamiltonian:[45]where is the effective electronic spin operator and I⃗ is the spin of nucleus M. The expectation value of this operator over pure spin functions, with spin quantization along the v-axis, gives the energy due the hyperfine interaction:This energy will be equal to the true hyperfine interaction energy,[46−48] obtained via a quantum mechanical description, EQM((I⃗). In other words, the result for the effective spin Hamiltonian can be mapped onto the results of the quantum mechanical Hamiltonian.[16] In order to determine an element of the hyperfine coupling tensor, the derivative with respect to the uth component of the nuclear spin is taken:

In the following, an appropriate quantum mechanical operator describing the hyperfine interaction will be derived starting from the relativistic Dirac Hamiltonian, with the electron–electron interaction given by the Coulomb operator:[49]where α⃗ and β are the Dirac matricesand σ⃗ is the vector consisting of the Pauli spin matrices:The nuclear potential in eq , Vnuc(i), is approximated by a finite nuclear charge distribution in the shape of a Gaussian function.[50]

To derive the operator for the hyperfine interaction, the magnetic field from the Mth nucleus is introduced in the Dirac Hamiltonian via the minimal coupling (using the cgs system of atomic units):[51]where A⃗ is the vector potential; within a point-like description of the magnetization distribution it is given bywhere μ⃗ is the magnetic moment of nucleus M given by μ⃗ = gμNI⃗, with g the nuclear g-factor and μN the nuclear magneton (μ = (2mc)−1).

Keeping only the term including A⃗ gives the 1-electron hyperfine interaction operator:and inserting the expression for the vector potential yields

In the case of variational wave functions (such as Hartree–Fock, DFT, CI, etc.) the derivative in eq can be translated into an expectation value using the Hellmann–Feynman theorem. In this work we employ the finite field method,[28] where the derivative is evaluated numerically. In the finite field method the perturbation operator is added to the zeroth-order Hamiltonian, eq , with a prefactor, λ, referred to as the field strength and proportional to I:

An element of the hyperfine coupling matrix can now be calculated as

The superscript (v) on the CC energy indicates the quantization axis of the total electronic angular momentum. This axis is in the present work controlled by taking advantage of the symmetry scheme employed by the Dirac program in which (for the symmetries considered here) the quantization axis is fixed along the z-axis.[52,53] ⟨S̃⟩ is simply the effective electronic spin and we will denote it S̃.

Due to the axial symmetry in diatomic molecules, the hyperfine interaction tensor can be described in terms of the parallel and the perpendicular components, denoted A∥ and A⊥. If the diatomic molecule is placed along the z-axis, A∥ and A⊥ can be calculated asand

Consequently, two separate coupled cluster calculations should be performed, the difference between them being the quantization axis of total electronic angular momentum. This scheme has been adopted from the approach presented in ref (54). In practice, the perpendicular component is obtained by placing the internuclear axis on either the x- or y-axis while the quantization axis of total electronic angular momentum is kept along the z-axis, effectively using the expression in eq . A similar scheme was recently presented in the framework of the complex generalized Hartree–Fock and Kohn–Sham methods.[55]

Computational Details

All the calculations were carried out with the DIRAC17 program package.[53] In addition to the relativistic four-component (4C) calculations also the exact two-component (X2C) method was employed.[56] The bond length of the BaF radical was taken from the NIST Chemistry WebBook and has the value of 2.162 Å.[57,58] For the two isotopes considered in this work, 133Cs and 137Ba, nuclear spins of 7/2 and 3/2 and magnetic moments of 2.582 μB and 0.937 μB, respectively, were taken from ref (59).

Basis Sets

We employ Dyall’s relativistic basis sets from the valence, vXz, and core–valence, cvXz, series, where X denotes the cardinal numbers double-, triple-, and quadruple-ζ.[60−62] The vXz basis sets include correlation functions (of up to d-, f-, and g-type for Cs and Ba) for the valence region which is defined as 5s5p6s6p. The cvXz basis sets include additional correlation functions (of up to f-, g- and h-type for Cs and Ba) for the core–valence region, which includes the 4d shell in addition to the 5s5p6s6p shells. The effect of adding particular types of tight functions, i.e., basis functions with large exponents, was investigated by adding functions in an even-tempered fashion.

Correlation Treatment

The unrestricted CC module (RELCC) of DIRAC was employed with different types of perturbative triples:[63] the widely used CCSD(T) method,[64] which includes some fifth-order triples contributions, the CCSD+T (also called CCSD[T]) method,[65] in which triples contributions only up to the fourth order are included, and the CCSD-T method,[66] where one further fifth-order triples diagram is added to the ones included in the CCSD(T) method.[63] The CCSD-T method is therefore formally the most complete method of the three, but its performance was shown to be very similar to CCSD(T).[32,66] In addition, we have employed the multireference Fock-space CC method (FSCC).[67,68] We have tested the (0,1) sector with varying size of the model space. In sector (0,1) a manifold of singly excited states are obtained by adding an electron to a closed shell singly ionized reference state. The additional electron can occupy those orbitals contained in the so-called model space. We will distinguish between two model spaces: A minimum model space (min) only including the valence orbital and an extended model space (ext) that includes the valence orbital as well as the five lowest virtual orbitals.

In both the single-reference CC and the FSCC calculations all electrons were included in the correlation calculation and consequently a high virtual space cutoff of 2000 au was used if not stated otherwise.

Finite Field Method

As a consequence of the introduction of the perturbation in eq , the total energy can be written as a Taylor series in λ:

The magnitude of λ should be chosen such that higher order terms will be negligible, i.e., E(λ) behaves linearly with small variations in λ. If indeed E(λ) is linear with respect to the variations in λ the two-point formula can be used to obtain the derivative:

With this two-point formula, any quadratic terms cancel out, resulting in an error proportional to λ2, as shown in ref (69) and the Supporting Information. Field strengths should be chosen large enough so that numerical instabilities are avoided and small enough so that higher order terms can safely be neglected. Therefore, a strict convergence criterion of 10–12 au for the CC amplitudes was used in the calculations.

Procedure

Since the HFS operator introduced above (eq ) is odd with respect to the time-reversal symmetry, it cannot be added directly on the DHF level, which in the DIRAC program is based on the Kramers-restricted formalism (krDHF). Instead, we add the operator on the CC level, which uses the unrestricted formalism. Consequently, both spin-polarization as well as correlation effects are accounted for by the CC iterations. In order to disentangle spin polarization and correlation effects, we also performed calculations on the Kramers-unrestricted DHF level (kuDHF) using the ReSpect program.[70] For a description of the kuDHF method, we refer to refs (71−73).

For clarity we outline the procedure of the calculation below. We note that the finite field scheme has long been available in the DIRAC program but has not, to our knowledge, been applied to HFS constants. In order to construct the HFS operator we simply employ operators from the catalogue of 1-electron operators included in the DIRAC program. The scheme is as follows:

Perform an unperturbed Kramers-restricted DHF calculation.

Carry out the integral transformation including integrals over the HFS operator, eq .

Determine the DHF energy in the presence of the field from the recomputed Fock-matrix. This will correspond to the Kramers-restricted DHF energy.

Perform two Kramers-unrestricted CC calculations in the presence of the positive and negative field to get the field dependent CC energies.

Calculate the numerical derivative of the CC energy using the two-point formula, eq .

Results and Discussion

Numerical Accuracy

Before turning to the effects of basis set, electron correlation, and relativity, we devote a section to the investigation of the numerical stability of the scheme presented above. In the case of the finite field method, special care must be taken to avoid numerical instabilities. For this purpose the X2C method and the vdz basis set have been used and only the parallel component, A∥, of the 137BaF HFS tensor has been considered, as the behavior is expected to be the same for the perpendicular component, A⊥.

In order to determine the appropriate field strengths to use with the finite field method, we investigated the dependence of the calculated HFS constants on the field strength. The HFS constants of 137BaF and 133Cs on the DHF, CCSD, and CCSD(T) level are shown in Table for the field strengths 10–9, 10–8, 10–7, 10–6, 10–5, 10–4, 10–3, 10–2, and 10–1 au. In all cases, the results for the lower field strengths of 10–9, 10–8, and 10–7 differ slightly from those obtained with the larger field strengths, indicating numerical instability. Whereas calculations with larger fields all yield the same values of the HFS constant (to the digits shown in the table) at the DHF level, the results on the CC level begin to deviate again at field strengths of ≥10–2. Note that the different dependence of the Hartree–Fock and CC results on the field strengths was also observed and discussed in detail in ref (74). The results for field strengths between 10–6 and 10–3 are stable for all methods, which indicates that the terms in the Taylor expansion (eq ) higher than quadratic are negligible (recalling the cancellation of quadratic terms by the 2-point formula). We have checked this by fitting the total energy as a function of λ to a third-order polynomial and found that the third order terms only become significant for field strengths above 10–3 au (see Supporting Information for further details). From the same fit the error due to neglecting the third order terms (by using the 2-point formula) at field strengths of 10–6 au can be estimated to be on the order of 10–10 au. We have thus chosen to use the two-point formula with a field strength of 10–6 au for all further calculations.

Table 1

Calculated A∥ and A Constants (MHz) of 137Ba in BaF and 133Cs for Different Field Strengthsa

	137BaF			133Cs
field	DHF	CCSD	CCSD(T)	DHF	CCSD	CCSD(T)
10–9	1650.2	2244.9	2244.9	1500.6	2114.8	2097.3
10–8	1644.3	2244.9	2230.0	1493.6	2110.4	2099.0
10–7	1645.3	2247.0	2233.3	1493.0	2109.6	2097.8
10–6	1645.2	2246.7	2233.2	1493.0	2109.5	2097.6
10–5	1645.2	2246.7	2233.2	1493.0	2109.5	2097.7
10–4	1645.2	2246.7	2233.2	1493.0	2109.5	2097.7
10–3	1645.2	2246.7	2233.2	1493.0	2109.5	2097.7
10–2	1645.2	2246.4	2232.9	1493.0	2109.2	2097.4
10–1	1645.2	2216.8	2203.1	1493.0	2087.3	2075.1

The calculations were performed using the X2C method and the vdz basis set.

It should be emphasized that the analysis described above should be performed for any new system in consideration. As an example take instead the 19F HFS constant in BaF, which is around 30 times smaller than the 137Ba and 133Cs HFS constants. The range of numerical instability is consequently larger (up to 10–6 au) for the Ba19F results and one would need to use larger field strengths (see Supporting Information).

To test the numerical accuracy further, we have performed a series of tests with the results listed in Table . The first test is related to the dependence of the CC HFS constants on the Hartree–Fock orbitals. We tested two different SCF convergence criteria of 5 × 10–9 and 1 × 10–8, resulting in a minor change of 0.05 and <0.01 MHz for BaF and Cs, respectively.

Table 2

Calculated A∥ and A Constants (MHz) of 137Ba in BaF and 133Cs for Various Computational Tests (Further Details in the Text)a

test		137BaF	133Cs
SCF convergence	1 × 10–8	2233.59	2098.10
	5 × 10–9	2233.54	2098.10
(SS|SS)	exclude	2233.54	2098.10
	include	2233.30	2097.86
screening	1 × 10–12	2233.54	2098.10
	1 × 10–15	2233.56	2098.10
	off	2233.56	2098.12

The calculations were performed using the X2C method and the vdz basis set.

Next we tested the effect of two computational approximations that are commonly employed to speed up the SCF calculations. The first is related to the inclusion of Coulomb integrals. The integrals involving only small-component wave functions, (SS|SS), have in all calculations been replaced by a simple Coulombic correction[75] and the effect of including them is here seen to be −0.24 MHz for both systems. This corresponds to less than 0.02% of the total values and is similar to that observed in previous studies of contact densities.[35] Second, we tested the effect of screening the 2-electron integrals used in the Fock matrix, that is, neglecting those estimated to be below a given threshold.[49] A threshold of 10–12 au is used as default in the DIRAC program and we find that turning the screening off (and thus including all 2-electron integrals) has a negligible effect of 0.02 MHz for both systems.

Using field strengths of 10–6 au and employing the approximations described above, we conclude that we can safely include four digits in the following discussions.

Basis Set Effects

Here we investigate the effect of the basis set on the calculated HFS constants. In order to reach the highest possible accuracy, we need to choose a basis set that is sufficiently converged with respect to additional functions. We consider the convergence sufficient when additional basis functions do not change the HFS constants by more than ∼0.5%, since we expect the total uncertainty of a few percent. At the same time the basis set should be small enough to allow for realistic CC calculations with large active spaces. The following basis set studies were carried out at the four-component CCSD level correlating all electrons and using a virtual cutoff of 2000 au, which will be justified in the section Correlation Effects.

In Table the HFS constants of 137Ba in BaF and 133Cs are shown with increasing quality of the valence and core–valence basis set series, vXz and cvXz (X = d (double), t (triple), q (quadruple)). For both series and both systems a converging behavior is observed upon increasing basis set quality, with the Cs results converging notably faster than the BaF results.

Table 3

Calculated A∥, A⊥, and A Constants (MHz) of 137Ba in BaF and 133Cs for Increasing Basis Set Qualitya

	137BaF				133Cs
	A∥	%(expb)	A⊥	%(expb)	A	%(expc)
vdz	2247	–5.4	2168	–5.8	2110	–8.2
vtz	2316	–2.5	2238	–2.7	2206	–4.0
vqz	2342	–1.4	2264	–1.6	2232	–2.9
s-aug-vqz	2342	–1.4	2265	–1.6	2232	–2.9
cvdz	2292	–3.5	2214	–3.8	2161	–6.0
cvtz	2363	–0.5	2285	–0.7	2264	–1.5
cvqz	2383	0.3	2305	0.2	2283	–0.7
aeqz	2386	0.4	2308	0.3	2287	–0.5
exp	2376(12)		2301(9)		2298.16

The calculations were performed using the 4C CCSD method. Deviation from the experimental values is also shown.

Reference (77).

Reference (78).

The addition of one diffuse function for each angular momentum to the vqz basis set, denoted s-aug-vqz, has negligible effect on the calculated HFS constants. This is as expected since the HFS constants describe the interaction of the unpaired electron with the Ba or the Cs nuclei and thus should not be strongly affected by the quality of the description of the region far away from the nuclei. Note that this is not the case for the HFS constants of excited states, where diffuse functions are of great importance.

The difference between the (c)vtz and (c)vqz results (of approximately 1%) indicates, however, that the basis set is not yet saturated with respect to this property. This can be attributed to the slow basis set convergence of the CC methods.[76] In contrast, previous studies using four-component DFT methods and the same basis sets showed convergence already at triple-ζ level for the HFS constants.[71,73]

In Table we also show the deviation of the calculated HFS constants from the experimental results.[77,78] For both systems the cvXz HFS constants are higher than the vXz ones, corresponding to a smaller deviation from experiment. On the quadruple-ζ level the difference between the vqz and the cvqz values is ∼2%. The cvXz basis sets include large exponent (tight) functions with high angular momenta, which are needed to correlate the 4d shell (in the case of Ba and Cs) which can be considered as the core–valence region. Since we are correlating all the electrons and considering a property that involves interaction between the valence electrons and the nucleus it is to be expected that core–valence correlation functions are needed for obtaining high accuracy results.

In Table we show the effect of adding tight functions of different symmetries individually to the vqz basis set (the corresponding exponents are listed in Table S4 in Supporting Information). Since the behavior of the parallel and perpendicular component of the 137BaF HFS tensor with respect to basis set is very similar we only considered A∥ in this case. The conclusion is that only the addition of tight f-functions has an influence on the calculated values. The addition of one tight f-function has the largest effect of 1.0% and 1.4% for 137BaF and 133Cs, respectively. The addition of another three tight f-functions has a smaller additional effect of 0.8% and 1.0% and further tight f-functions are not expected to change the results by more than 0.2%.

Table 4

Calculated A∥ and A Constants [MHz] of 137Ba in BaF and 133Cs with Different Tight Functions Added to the vqz Basis (Corresponding Exponents in Table S4)a

	137BaF		133Cs
X	A∥		A
vqz	2342	0.0	2232	0.0
+s	2342	0.0	2231	0.0
+p	2342	0.0	2232	0.0
+d	2342	0.0	2232	0.0
+f	2366	1.0	2262	1.4
+2f	2376	1.4	2274	1.9
+3f	2380	1.6	2281	2.2
+4f	2383	1.8	2285	2.4
+g	2343	0.0	2232	0.0
+h	2343	0.0	2232	0.0

The calculations were performed using the 4C CCSD method. The effect (in %) with respect to the vqz basis is also shown.

As the cvqz basis set differs from the vqz basis set by three tight f-, two tight g-, and one tight h-functions, we can conclude that the differences between the vqz and cvqz results are governed by the addition of the tight f-functions. To test that the cvqz is indeed converged with respect to the addition of tight functions, we used the all-electron quadruple-ζ basis (aeqz) set which includes correlation functions for all shells, resulting in a minor increase in the HFS constant of ∼0.2%. If not stated otherwise, we have thus chosen to use the cvqz basis sets in our further investigations of other computational parameters.

It has been shown previously that the addition of tight s-functions to standard correlation consistent basis sets is necessary to accurately calculate the HFS constants.[79,80] This is not the case here as seen in Table , indicating that the size of the Dyall vqz basis set in the vicinity of the nucleus is already sufficient.

Correlation Effects

Table and Figure contain the HFS constants of 137BaF and 133Cs, obtained at different levels of theory. In addition to the total HFS constants, the correlation contribution compared to the krDHF result is shown explicitly along with the deviation from experiment.

Table 5

Calculated A∥, A⊥, and A Constants (MHz) of 137Ba in BaF and 133Cs at Different Levels of Correlationa

137BaF
	A∥	Δ	%(expb)	A⊥	Δ	%(expb)
krDHF	1598	0	–32.8	1553	0	–32.5
kuDHFd	1905	307	–19.8	1817	260	–21.0
CCSD	2383	785	0.28	2305	752	0.19
FSCCSD min	2399	801	0.96	2323	770	0.94
FSCCSD ext	2403	806	1.16	2328	775	1.16
CCSD+T	2425	827	2.06	2350	797	2.14
CCSD(T)	2358	760	–0.77	2282	729	–0.85
CCSD-T	2365	767	–0.45	2288	735	–0.56

The cvqz basis sets were used in the calculations.

Reference (77).

Reference (78).

Results obtained with the ReSpect program.[70,71]

Figure 1

Calculated A∥, A⊥, and A constants (MHz) of 137Ba in BaF and 133Cs at different correlation levels, compared to experiment. The shaded areas indicate experimental uncertainties.

As expected, the lack of correlation treatment as well as of spin polarization in the krDHF method results in a significant underestimation of more than 30% compared to the experimental results. The inclusion of spin polarization in the kuDHF method leads to a significant increase in the HFS constants resulting in a deviation around 20%. However, one certainly needs to go to the CC methods for high accuracy.

With the CCSD method the HFS constants are thus significantly higher, resulting in a deviation from experiment of less than 1%. The multireference Fock-space CC method (FSCCSD) produces results between the CCSD and CCSD+T values, which is due to the fact that the FSCCSD method takes part of higher order contributions (beyond the double excitations of CCSD) into account due to its multireference formalism. Extending the model space used with FSCCSD (FSCCSD ext, see the section Correlation Treatment for a description of the employed model spaces) has a negligible effect, indicating good description of the two systems by a single reference determinant, 2Σ1/2 in the case of BaF and 2S1/2 in the case of Cs.

The inclusion of perturbative triples has a small effect, with the CCSD+T results slightly overestimating and the CCSD(T) and CCSD-T slightly underestimating the experimental values (see inset of Figure ). A similar nonsystematic behavior was observed in ref (32) for electric field gradients. However, the present findings are unusual in that the fluctuations in the size of the perturbative triples contributions obtained with the different approximations are comparable with their total values (that is, the difference between the CCSD and CCSD+T/(T)/-T results). For the effective field gradients[32] and in the recent studies of various P- and P,T-odd interaction constants,[37,38] these fluctuations were significantly smaller than the total contribution of the perturbative triple excitations.

Our results indicate that the triple excitations are more important for the HFS constants than for the other properties mentioned above. This has been recognized in the past, by, for example, Safronova et al.,[8] or more recently by Tang et al.,[81] who identified this issue from the relatively large difference between the linearized and the full CCSD method. Consequently, we choose to continue our analysis with CCSD and to base our recommended values and uncertainty estimates on this method.

The correlation contributions to the HFS constants are almost identical for A∥ in BaF and A in Cs whereas the correlation contribution to A⊥ in BaF is slightly lower. It is interesting to note that the trends and differences between the different methods are very similar in BaF and Cs, Figure . This indicates that the two system have a similar electronic structure. In BaF one of the two valence electrons of Ba is participating in the bonding to F leaving a Ba + like system, which is iso-electronic to the Cs atom.

The results presented until now have included correlation of all the electrons and a cutoff of 2000 au of the virtual correlation orbitals. As shown for example in ref (82), a high virtual cutoff is needed in order to capture the correlation contributions to HFS constants associated with the core electrons. In Figure we present in detail the dependence of the HFS constants on the virtual space cutoff when correlating all electrons in BaF and Cs. In both cases only specific virtual orbitals have a significant influence on the correlation contribution to the HFS constants. Inspection of the orbitals in question (see Supporting Information) reveals that the contributing orbitals are all of s-function character (s-functions of Ba in the case of BaF). From the deviation with respect to results obtained when all the virtual orbitals were included in the correlation space (designated “no cutoff” on the Figure y-axis) it can be seen that choosing a cutoff of 2000 au will result in an underestimation of the HFS constants of approximately 0.5%. Since this uncertainty is smaller than the expected uncertainty of the method we choose to proceed with a cutoff of 2000 au

Figure 2

Calculated A∥ and A constants (MHz) of 137Ba in BaF and 133Cs at the CCSD/vtz level for different virtual space cutoffs. See text for further details.

Relativistic Effects

So far we have presented results on the four-component Dirac–Coulomb (DC) level of theory. The last part of this analysis is dedicated to the investigation of the dependence of the calculated HFS constants on the treatment of relativistic (and related) effects. The results obtained using different models are shown in Table .

Table 6

Calculated A∥, A⊥, and A Constants (MHz) of 137Ba in BaF and 133Cs at Different Levels of Treatment of Relativistic Effectsa

	137BaF		133Cs
	A∥	A⊥	A
CCSD DC	2383	2305	2283
CCSD X2C	2382	2305	2283
CCSD DCG	2382	2305	2282
CCSD PN	2414	2337	2312

The cvqz basis sets were used in the calculations.

As expected, the X2C and DC Hamiltonians give practically identical results, confirming the excellent performance of the former.

In the DC Hamiltonian the 2-electron interaction is approximated by the Coulomb potential, which can be considered as a nonrelativistic description (it is instantaneous and not Lorentz invariant). For a proper relativistic description of this interaction one needs to turn to the theory of quantum electrodynamics (QED), where one takes into account the finite speed of light resulting in a noninstantaneous interaction. The lowest order one-photon exchange interaction in the static approximation can be derived in the Feynman gauge or the Coulomb gauge, referred to as the Gaunt and Breit interactions, respectively.[83] Whereas the Breit interaction is correct to , the Gaunt interaction is correct to and simpler to implement and calculate. The current implementation allows us to include the Gaunt interaction on the DHF level (DCG); these results are shown in Table . We observe a negligible effect of the Gaunt contribution of ∼−1 MHz on the HFS constants. Previous studies on 133Cs have considered the Gaunt[9] or the full Breit interaction[8,10,84,85] at different stages of the calculations. For a thorough comparison and discussion of some of these efforts we refer to ref (10). In comparison to the majority of the results (4.87 MHz,[10] 5.0 MHz,[9] and 6.00 MHz[85]), we, however, predict the wrong sign as well as too small an effect for the Gaunt interaction contribution, which might be due to several factors: first of all, we calculate the Gaunt contribution on the DHF level only, lacking any Gaunt contribution on the correlated level. Second, we employ the restricted DHF formalism, which might lack relaxation effects. Indeed, the negative Breit contribution obtained in ref (8) was attributed to the neglect of relaxation effects due to the perturbative approach.

Finally, we test the dependence of the HFS constants on the employed nuclear model. In Table we present results obtained using a point-like description of the nuclear charge (PN). Despite the seemingly big physical difference between the point-like and Gaussian description of the nuclear charge, the effect on the calculated HFS constants is relatively small (1.3% for A∥ in BaF and A in Cs and 1.4% for A⊥ in BaF). Nonetheless, the Gaussian model should be employed if high accuracy is desired. In previous studies on the DFT level,[71,86] the effect of the finite size of the nuclear charge distribution was found to be ∼1% for Zn HFS constants, ∼1.5% for Cd HFS constants and as large as ∼10–15% for Hg HFS constants.

The authors of ref (71) also investigated the effect of a Gaussian description of the nuclear magnetic moment distribution, which turned out to be negligible for lighter elements and as large as ∼2% for Hg. This effect was also studied by Ginges et al.[85] who found contributions ranging from 0.18(15)% for 133Cs to 4.35(131)% for 225Ra, which shows that a finite distribution of the magnetic moment should be included if a small uncertainty is desired for the HFS constants of the sixth row elements. The fact that we neglect this effect in the present calculations is one of the main sources of uncertainty, especially for 137BaF (see the section Uncertainty Estimation).

Uncertainty Estimation

On the basis of the investigations presented in the previous sections, we consider the results on the CCSD DC/cvqz level to be our recommended values. On this level of theory the convergence with respect to basis set was sufficient and the correlation treatment was the most reliable.

In addition to the comparison with experimental results, we perform an uncertainty analysis based purely on theoretical considerations. In cases where no experimental data are available, a theoretical uncertainty estimate is crucial for direct applications of the calculated properties in experimental research. Here we follow a similar procedure to that in our previous work on symmetry breaking properties.[37,38] In this scheme we estimate the error that is introduced by the different approximations employed in the treatment of the basis sets, electron correlation, relativistic effects and nuclear description. The individual sources are added in a quadratic fashion, which assumes the considered effects to be independent. While it is well-known that computational effects, associated with for example the treatment of electron correlation and choice of basis set, depend strongly on each other, this assumption can be considered good enough as long as the individual uncertainties are small which is the case here. These sources of uncertainty are presented in Table and discussed in the following.

Table 7

Summary of the Sources of Uncertainty (MHz) of the Calculated A∥, A⊥, and A Constants (MHz) of 137Ba in BaF and 133Cs

source	137BaF		133Cs
	δA∥	δA⊥	δA
basis set
quality	20.00	20.00	19.0
tight functions	3.00	3.00	4.00
diffuse functions	0.00	1.00	0.00
correlation
higher order	–120.00	–124.00	–120.00
virtual cutoff	8.18	8.18c	12.78
relativistic effects
Breit	5.72a	5.53a	6.00b
QEDVP+SE	–10.01a	–9.68a	–10.30b
Bohr–Weisskopf	–39.56a	–38.26a	–7.60b
quadratic sum	128.74	132.07	123.05
%	5.40	5.73	5.28

Based on 135Ba+ results from ref (85).

Taken directly from ref (85).

Used A∥ results.

Basis Set

In section Basis Set Effects we investigated the effect on the HFS constants of increasing the basis set size in three aspects; the addition of tight functions, diffuse functions and the general quality. We finally chose to use the cvqz basis set and we estimate the uncertainty that is introduced by truncation at the quadruple-ζ level to be not larger than the difference between the cvtz and cvqz results. The effect of adding additional tight (aeqz) and diffuse (s-aug-vqz) functions turned out to be very small but we include them here for the sake of completeness. Adding all three effects together amounts to 23, 24, and 23 MHz for both A∥ and A⊥ in 137BaF and A in 133Cs, which corresponds to a bit more that 1%.

Electron Correlation

In our previous studies we used the spread in the perturbative triples results (i.e., the difference between the CCSD+T and CCSD-T results) times 2 as an estimate for the order of magnitude of the missing higher order correlation contributions.[37,38] In both cases this was close to half of the difference between CCSD and CCSD(T). However, in the case of the HFS constants the difference between CCSD+T and CCSD-T is ∼60 MHz for both systems, about 3 times larger than the difference between CCSD and CCSD(T). This is an indication that higher order correlation contributions are more important in the case of HFS constant. As a conservative estimate, we use again the spread in the perturbative triples results multiplied by 2, which is the major source of uncertainty and contributes ∼5% in both cases.

In the section Correlation Effects we found that neglecting the virtual orbitals above 2000 au introduces an error of ∼0.5% and we add this contribution to the uncertainty estimate.

Relativistic Effects (Breit and QEDVP+SE)

In order to estimate the magnitude of the higher order relativistic corrections to the 2-electron interaction, we rely on previous works and in particular on the recent study by Ginges et al.[85] who systematically investigated various contributions to the ground state HFS constants of a few atoms and ions.

A thorough discussion on the previous calculations of the Breit contribution to the HFS constant in 133Cs can be found in ref (10) where also the, at the time, most rigorous calculation of the Breit contribution at the level of third-order many-body perturbation theory (MBPT) was presented being 4.9 MHz. In the recent study by Ginges et al.[85] this contribution was estimated to be 6.0 MHz at the level of the random phase approximation (RPA). We use the larger value of Ginges et al. to estimate the effect of neglecting the Breit interaction.

To our knowledge, no study of the Breit contribution to the 137BaF HFS constant was published to date. Due to the similar electronic structure and nuclear charge of 137BaF and 133Cs the Breit contribution is expected to be similar and we could use the 133Cs results as an estimate for the effect in 137BaF. We choose instead to estimate this effect from the result in ref (85) for the 135Ba+ HFS constant. The electronic structure of the Ba+ ion is a good approximation to that in BaF, where one of the two valence electrons of Ba participates in the bonding to F leaving Ba effectively with a positive charge. The isotope effect on the Breit contribution is negligible. The Breit contribution was determined in ref (85) to be 0.24% of the total HFS constant of 135Ba+. Taking this to be representative for the 137BaF HFS constant, we estimate the Breit contribution as 5.72 MHz for A∥ and 5.53 MHz for A⊥, which indeed is very similar to that in 133Cs.

For higher order corrections to the 2-electron interaction, one has to turn to quantum electrodynamics (QED) where the lowest order diagrams (beyond Breit) are the single photon one-loop diagrams, namely, the vacuum polarization and the self-energy, QEDVP+SE.

Two predictions of the QEDVP+SE contributions to the HFS constant of Cs are available. One is from Sapirstein et al.[87] of −9.7 MHz, and the other from Ginges et al.[85] of −8.8(15) MHz, which agree within the uncertainty provided for the latter. As an estimate, we choose the latter value, including the provided uncertainty. For 135Ba+ Ginges et al. predicted −0.38(4)%, which translates to −10.01 and −9.68 MHz for A∥ and A⊥ in 137BaF.

Bohr–Weisskopf Effect

Finally, we consider the Bohr–Weisskopf effect, which accounts for the finite distribution of the nuclear magnetization compared to a pointlike model employed in this work. Again we use the results from ref (85), which, unlike the Breit and the QEDVP+SE effects, turn out to be quite different for the two systems, i.e., −0.18(14)% for 133Cs and −1.26(38)% for 135Ba+. This difference originates from the different nuclear properties of the two isotopes. We have also investigated this effect on the kuDHF level using the approach described in ref (71), where the nuclear magnetic moment distribution is described by a simple Gaussian function. The results are shown in Supporting Information and the effect is ∼−0.8% for both systems. While the sign and order of magnitude are similar to the higher level results from ref (85), this method is not able to take the different nuclear structures into account. However, this method will be very useful in cases where higher level reference results are not available and in particular for larger molecules.

The similar nuclear properties of the 135Ba and 137Ba isotopes result in a similar Bohr–Weisskopf effect[88] and we use the estimate for 135Ba+ also from ref (85) in our uncertainty estimate. We note that besides nuclear structure the Bohr–Weisskopf effect also strongly depends on the electronic state of the system, which was recently demonstrated by Prosnyak et al. for Tl.[89]

Comparison with Previous Studies

Before we conclude, we compare our results with earlier theoretical values and with experimental results. Since the Gaunt contribution was seemingly unreliable, i.e., predicting the wrong sign, and the perturbative triples contributions seemed unreliable due to their relatively large spread, we choose the DC CCSD results (using the cvqz basis set) to be our best estimate for the HFS constant in these two systems, with the associated uncertainties presented in Table .

For both systems the deviation of the DC CCSD results from the experimental values is below 1%, as can be seen in Tables and 9. This deviation is well below the estimated uncertainty of >5%, Table . It illustrates the conservative nature of our error estimate, in particular in the higher order correlation corrections, but it is also a result of cancellation between the uncertainties stemming from basis set, correlation, and Breit interaction (positive) and the QEDVP+SE and Bohr–Weisskopf effects (negative).

Table 8

A∥ and A⊥ of 137Ba in BaF (MHz)

	137BaF
method	A∥	%(exp)	A⊥	%(exp)
GRECP SCF-EO[90]	2264	–4.71	2186	–5.00
GRECP RASSCF-EO[90]	2272	–4.38	2200	–4.39
DF RASCI[91]	2240	–5.72	2144	–6.82
DF MBPT[91]	2314	–2.61	2254	–2.04
DC CCSD (this work)	2383(129)	0.29	2305(132)	0.17
exp[77]	2376(12)		2301(9)

Table 9

A of Cs in MHza

method	133Cs	%(exp)
MBPTb+B[7]	2291.00	–0.31
SDpT+B[8]	2278.5	–0.85
MBPTb[7]+B[10]	2295.87	–0.10
MBPTb+OE+G[9]	2302	0.17
CCSDvT[11]+B[10]+QEDVP+SE [87]	2306.6	0.36
CCSD (ECC)[26]	2179.1	–5.18
CCSD (Z-vector)[27]	2218.4	–3.47
MBPTb+B+QEDVP+SE [85]	2294.4	–0.16
CCSD (LCCSD)[81]	2345.9	2.08
CCSD (finite field, this work)	2283(123)	–0.66
exp[78]	2298.16

All methods employed the four-component formalism. +B and +G denote the inclusion of the Breit and Gaunt interaction, respectively. For the CCSD methods the procedure used to extract the HFS constant is given in parentheses.

MBPT has been used as a general term for atomic many-body methods. While the MBPT results were all obtained using Brueckner orbitals in the evaluation of HFS matrix elements (at the RPA level) there are some smaller differences between the methods.

BaF

Two previous studies have reported calculations of the 137BaF HFS constant; these results are presented in Table . The first study by Kozlov et al.[90] reported results obtained with the self-consistent field (SCF) and restricted active space SCF (RASSCF) methods with and without core-polarization included with the aid of effective operators (EO). The effect of including core polarization (∼780 MHz for A∥ and ∼740 MHz for A⊥) was seen to be very similar to the effect of going from SCF to CCSD discussed in section Correlation Effects. Furthermore, the RASSCF-EO show little difference to SCF-EO, which agrees with the small difference between CCSD and FSCCSD. The restricted active space configuration interaction (RASCI) result of Nayak et al.[91] is very similar to the (RAS)SCF-EO results, both underestimating the HFS constant by about 5% compared to the experimental value. The use of MBPT offers a significant improvement compared to the RASCI results.

From the results listed in Table the present DC CCSD result has the smallest deviation from the experimental value and offers an improvement of accuracy compared to the earlier investigations.

Cs

The HFS constant of Cs has been studied extensively due to its relevance for atomic parity violation experiments.[5,6] Interpretation of such experiments requires sub 1% accuracy for the theoretical predictions. As can be seen from Table , this goal has been achieved by several groups over the years using various many-body methods.[7−11,85] Most of the results with less than 1% deviation from experiment were obtained with atomic codes, where use of the radial symmetry can practically eliminate basis set errors. Another feature of these results is that they all include a subset of triple excitations as well as estimates for the Breit and/or QEDVP+SE corrections. Therefore, while the present DC CCSD values have a similar error with respect to experiment, a direct comparison with the earlier high accuracy studies is not meaningful.

In recent years Sasmal and co-workers have reported the HFS constants of a large set of atoms and molecules on the CCSD level using the extended CC (ECC) and Z-vector frameworks.[26,27] The ECC uses a variational CC ansatz that allows for calculating HFS constants as expectation values. The Z-vector technique, however, is a way to evaluate the energy derivative of nonvariational CC energies. Due to the cumbersome truncation scheme in the case of ECC the Z-vector approach is expected to perform better. Indeed, the deviation with respect to experiment is smaller for the Z-vector result compared to the ECC result but still significantly larger than the aforementioned many-body methods. There can be several reasons for this; first of all, these results were obtained with molecular codes that would suffer from similar basis set uncertainties as presented in this work. Second, the ECC as well as the Z-vector results were obtained with a virtual cutoff of 60 and 40 au, respectively. This cutoff corresponds to the first few points in Figure , which indeed leads to an underestimation of ∼3%. The advantage of the present finite field approach over the ECC and Z-vector methods is that it allows for the inclusion of perturbational triples. Even though the perturbational treatment of the triple excitations turned out to be too unstable to include in the recommended value, it still provides an important contribution to the uncertainty estimation.

Recently, an additional study on the DC CCSD level was reported by Tang et al.[81] In their approach the linearized expression for the CCSD expectation value was employed while the amplitudes were obtained from a CCSD calculation taking all terms into account. The overestimation of ∼2% was attributed to the missing nonlinear terms in the expectation value expression.

Conclusion

We calculated the HFS constants of 137BaF and 133Cs on the relativistic coupled cluster level using the finite-field method as a straightforward way to evaluate the energy derivative. This scheme has been previously applied to various properties, but the present work is the first application to HFS constants. Consequently, a detailed investigation of computational parameters has been performed and presented. The effect of including different types of perturbative triples on the calculated HFS constants was seen to be more irregular than in the previous studies. We thus expect triple excitations to be important and conclude that a perturbational treatment is insufficient.

On the basis of the computational investigations, a transparent theoretical uncertainty estimate has been performed. Because of the irregular behavior of the perturbative triples, the largest contribution to the uncertainty estimate comes from the higher order correlations. Higher order relativistic as well as nuclear magnetization distribution effects were included in the estimate by using results from the literature. The estimated uncertainties amounted to 129 MHz (5.4%) and 132 MHz (5.7%) for A∥ and A⊥ in 137BaF and 123 MHz (5.28%) for 133Cs. These uncertainties are notably larger than those predicted for the P,T-odd interaction constants (∼2%) that were obtained using the same scheme as in the present work.[37,38,92]

The estimated uncertainties were found to be well above the deviation from experimental results which for both systems was below 1%. This discrepancy is partly due to the conservative nature of the uncertainty estimate (especially in the case of the higher order correlation effects), but it also reflects a fortunate cancellation of the missing contributions. An important task for the future is consequently to improve the description of higher order correlations that would enable more reliable uncertainty estimates.