Molecular enhancement factors for P, T-violating eEDM in BaCH$_3$ and YbCH$_3$ symmetric top molecules

High-precision tests of fundamental symmetries are looking for the parity- (P), time-reversal- (T) violating electric dipole moment of the electron (eEDM) as proof of physics beyond the Standard Model. Particularly, in polyatomic molecules, the complex vibrational and rotational structure gives the possibility to reach high enhancement of the P, T-odd effects in moderate electric fields. Additionally, it is possible to increase the statistical sensitivity by using laser cooling. In this work, we calculate the P, T-odd electronic structure parameters $W_\mathrm{d}$ and $W_\mathrm{s}$ for the promising candidates BaCH$_3$ and YbCH$_3$ for the interpretation of future experiments. We employ high-accuracy relativistic coupled cluster methods and systematically evaluate the uncertainties of our computational approach. Compared to other Ba- and Yb-containing molecules, BaCH$_3$ and YbCH$_3$ exhibit larger $W_\mathrm{d}$ and $W_\mathrm{s}$ associated to increased covalent character of the M--C bond. The calculated values are $3.22\pm 0.11 \times 10^{24}\frac{h\text{Hz}}{e\text{cm}}$ and $13.80\pm 0.35 \times 10^{24}\frac{h\text{Hz}}{e\text{cm}}$ for $W_\mathrm{d}$, and $8.42\pm0.29$~$h$kHz and $45.35\pm1.15$~$h$kHz for $W_\mathrm{s}$, in BaCH$_3$ and YbCH$_3$, respectively. The robust, accurate, and cost-effective computational scheme reported in this work makes our results suitable for extracting the relevant fundamental properties from future measurements and also can be used to explore other polyatomic molecules sensitive to various violations of fundamental symmetries.


I. INTRODUCTION
The Standard Model of particle physics (SM) is known to be the most successful theory in describing the universe at the smallest scale. This model predicts all the known fundamental particles and explains their interactions via three of the four fundamental forces. However, despite its successful descriptions, SM is known to be an incomplete theory. Several well established experimentally observed facts, such as the matter-antimatter asymmetry and the existence of dark matter [1] and dark energy [2] are not described by the SM. In particular, the matter-antimatter asymmetry requires an amount of charge-parity (CP) violation incompatible with the SM [3]. The incompleteness of the SM is both an incentive and an opportunity to look for physics beyond the Standard Model (BSM), also known as new physics.
High-precision tests of fundamental symmetries are a very effective means of probing BSM physics [4]. Specifically, precision experiments in the sub-microhertz level using atoms and molecules are searching for the electric * lukas.f.pasteka@uniba.sk dipole moment of the electron (eEDM). The eEDM violates both time-reversal symmetry (T) and parity symmetry (P) and, assuming the CPT theorem, the eEDM thus also violates CP symmetry. In the SM, the eEDM is highly suppressed, and the predicted value is far too small to be measured using current experimental techniques. This has been estimated to be in the order of magnitude of |d SM e | < 10 −40 [5]. On the other hand, the BSM theories predict values in the experimental reach [6,7], and a measurement of a non-zero value would be an incontrovertible proof of new physics [8].
Presence of the eEDM induces an EDM on paramagetic molecules [9,10] which is enhanced due to the internal electric fields. It has been previously shown that this enhancement grows with the atomic number as Z 3 [11] making systems containing heavy atoms ideal for measuring the eEDM. Additionally, the experimental signal is also enhanced when using close-lying opposite parity states -more information can be found in the work of Sandars [12,13]. In atoms, opposite-parity electronic states are split by ∼ 2 eV, while, in molecules, opposite-parity rotational states are split by ∼ 10 −5 eV, typically. Consequently, some molecular states can offer a dramatically higher enhancement than the enhance-ment found in atoms [14,15]. Motivated by such enhancements, numerous experiments are being developed using molecules containing heavy atoms. The first experiments using diatomic molecules were performed on TlF in Oxford [14,[16][17][18], and the current most stringent result has been set in ThO by ACME collaboration; |d ThO e | < 1.1 × 10 −29 e·cm [19]. Other investigated diatomic molecules include YbF [20], HfF + [21], ThO [19], RaF [22,23], and BaF [24].
Some diatomic molecules, such as PbO [25], ThO [19], and HfF + [21], have almost degenerate opposite parity (excited) eigenstates, called Ω-doublets, which are used to measure the eEDM. These parity doublets have a small splitting which makes it possible to fully mix (or polarize) them in moderate electric fields. Additionally, the use of Ω-doublets give the possibility to cancel many systematic effects. Unfortunately, due to the complex structure of the Ω-doublets states, it is not possible to take advantage of laser cooling to improve the experimental sensitivity. On the other hand, measurements of the eEDM in diatomic molecules in their ground state, such as BaF [24], YbF [20,26], and RaF [23], cannot make use of the Ωdoublet states.
Polyatomic molecules emerge as good candidates for eEDM experiments [27,28]. In contrast to diatomic molecules, in polyatomic molecules it is possible to measure the eEDM in the long-lived close-lying opposite parity eigenstates (K-doublets). In addition, previous works [29][30][31][32][33][34] have discussed the feasibility of laser cooling of metal-containing molecules composed of a single metal atom bound to a single molecular ligand, such as the symmetric top molecules BaCH 3 and YbCH 3 . This is also supported by that fact the lighter CaCH 3 and MgCH 3 exhibit quasi-diagonal Frank-Condon factors [31].
The enhancement of the eEDM interaction in atoms and molecules is usually expressed in terms of the electronic structure parameter W d . The scalar-pseudoscalar (S-PS) interaction between the electrons and the nucleons is also a P,T-violating interaction expressed in terms of the parameter W s [35] and is similarly enhanced in atoms and molecules. Both parameters are needed to connect the measured energy shift due to P,T-violation to the eEDM value. Since W d and W s can not be measured experimentally, it is necessary to use electronic structure methods to calculate them.
In this work, we report the W d and W s parameters for the promising eEDM candidates for future experiments, the BaCH 3 and YbCH 3 molecules. Since the enhancement factor depends mainly on the identity of the heavy atom, it is expected that BaCH 3 and YbCH 3 have similar enhancement factors comparable to other isolectronic molecules relevant for eEDM experiments, such as BaF, BaOH, and YbF, YbOH [36][37][38][39], respectively. However, unlike in the diatomic molecules BaF and YbF, the symmetric top molecules BaCH 3 and YbCH 3 have long-lived K-doublets accessible to experimental measurement with an even smaller splitting than in the linear polyatomic BaOH and YbOH molecules, typically ≤ 1 MHz [27].
The use of these K-doublets to measure the eEDM in BaCH 3 and YbCH 3 makes it possible to access large enhancement of the eEDM using moderate electric fields, reach a high experimental sensitivity, and avoid many systematic effects.
To calculate the W d and W s parameters, we employ high accuracy single reference and Fock-space coupled cluster methods and we explore the effect of the different computational factors on the calculated values. The employed methodology allows us to estimate the uncertainty in the calculated values of W d and W s . The accurate and robust computational scheme established in this work may be extended to other polyatomic molecules sensitive to parity-violating effects.

II. METHODOLOGY
The eEDM operator, H eEDM , can be written in terms of a one-body operator [40] (here, and throughout the rest of this section, atomic units are used), where γ 0 , γ 1 , γ 2 , and γ 3 represent the Dirac matrices, γ 5 = iγ 0 γ 1 γ 2 γ 3 , p i is the momentum of the electron i, and c is the speed of light. The S-PS interaction can be expressed as where G F is the Fermi constant (2.2225 × 10 −14 a.u.), Z N is the atomic number of the nucleus N , and ρ(r iN ) is the nuclear charge distribution. In equations (1) and (2), d e and k s parameterize the eEDM and S-PS interaction, respectively.
The total Hamiltonian H is expressed as a sum of a zeroth-order Hamiltonian H (0) and a perturbation H k , regulated by the field strength parameter λ k , In our case, H (0) is the unperturbed molecular Dirac-Coulomb Hamiltonian, where α i and β i are the Dirac matrices and V nuc is the Coulomb potential. Considering the perturbations H eEDM and H S-PS , and applying the Hellman-Feynman theorem, the W d and W s coupling constants can be calculated as the first derivatives of the energy with respect to λ k , (6) and where Ψ (0) Ω is the unperturbed ground state electronic molecular wave function and Ω is the projection of the total electronic angular momentum on the internuclear axis.
The W d and W s constants are combined in the P,T-odd effective Hamiltonian, H P,T-odd [35], where S is the effective spin and n is a unit vector oriented along the internuclear axis. Therefore, a measurement of the P,T-violating energy difference on a single molecule provides us with the combined d e and k s and the disentanglement of the two effects requires experiments on different molecules [47].

III. RESULTS
In this work we assay a cost-accuracy balanced methodology that allows us to study polyatomic molecules promising for precision experiments, such as BaCH 3 and YbCH 3 , at the high accuracy coupled cluster level. We study the effect of the treatment of relativity and electron correlation and the choice of the basis set on the optimised molecular geometries (section III A) and on the calculated W d and W s parameters (section III B). Finally, we estimate the uncertainty of the predicted W d and W s based on an extensive computational study within the presented methodology (section III C).
All the calculations were carried out using a modified version of the Dirac 2019 program [48,49], except for the scalar relativistic (SR) calculations, where the CFOUR program [50,51] was employed. If not stated otherwise, the default settings of the corresponding codes were used. We applied both the single reference coupled cluster with single, double and perturbative triple excitations (CCSD(T)) [52] and the Fock-space coupled cluster (FSCC) with single and double excitations approach using sector (0,1) [53]. We used the uncontracted Dyall's relativistic basis sets [54][55][56] and the contracted atomic natural orbital correlation consistent ANO-RCC basis sets [57][58][59][60] of double-, triple-, quadruple-, and quintuplezeta (for the ANO-RCC basis sets only) cardinality.

A. Geometry optimisation
The spectroscopic properties and the geometries of the polyatomic molecules that are considered to be good candidates for measurements of P,T-violating phenomena are usually not known a priori, meaning that any computational study of these systems should begin with a geometry optimisation.
Due to the large number of electrons, geometry optimisation of polyatomic molecules at the 4-component coupled cluster level requires impractically high computational resources. We thus look for a compromise that allows sufficiently accurate calculations at a realistic computational cost. In the Appendix VI A we describe the effects of the relevant computational parameters on the optimised geometry and use these investigations to select the optimal computational approach. After evaluating the effect of treatment of relativity, use of contracted vs. uncontracted basis set, and electron correlation, we conclude that a suitable methodology for reliable geometry optimisations of small polyatomic molecules is the combination of the SR-CCSD/CCSD(T) approach with the contracted ANO-RCC basis sets and correlating the valence and core-valence ((n − 1) and (n − 2)) electrons.
Use of large basis sets is indispensable for taking full advantage of high accuracy correlation methods such as coupled cluster and for getting a good quality description of the system. Table I and Figure 1 show the convergence of the obtained geometry of the two molecules with the increase of the cardinality of the ANO-RCC basis set. The optimised geometry of BaCH 3 (converged at the 5z basis set level) is close (0.3 -0.7% difference) to the experimental geometry estimated from rotational constants in a gas-phase detection experiment (Ba-C bond distance 2.557 -2.570Å, BaCH angle 103 • -108 • , both assuming C-H bond distance 1.1Å) [61]. This slight underestimation is most likely due to the combined effects of scalar relativity and basis set contraction, as can be seen from Table X in Appendix VI A. We evaluate the effect of this discrepancy on the W d and W s parameters and include it as a source of uncertainty in our reported values (see Section III C).
In YbCH 3 , we found a faster convergence than in BaCH 3 , and thus used the geometry obtained with the QZ basis set. The optimised geometries used for the calculations of the W d and W s parameters in the following sections are presented in bold font in Table I.

B. Enhancement factors: Computational parameters
Both the W d and W s factors are purely relativistic properties and, according to Schiff's theorem, in the nonrelativistic regime the atom or the molecule would not have an EDM even if the electron did [62]. However, Schiff's theorem is not valid in the relativistic regime, and atoms and molecules may express a non-zero EDM  [12]. Therefore, calculations of the W d and W s factors should be carried out in a relativistic framework.
In this work, we use the relativistic Dirac-Coulomb 4-component Hamiltonian combined with the singlereference coupled cluster approximation [52] for the calculation of the W d and W s factors in BaCH 3 , and the multireference Fock-space coupled cluster approach [53], for the corresponding calculations in YbCH 3 .
Earlier calculations of the electric-field gradients in YbF [63] and of W d in YbF [37] and YbOH [38] have shown that Yb-containing molecules may present a multireference character in their ground states. In this work, we found, through the T1 diagnostic [64], that the ground state of YbCH 3 does, indeed, could indeed benefit from Fock space coupled cluster method. In the FSCC(0,1) approach, the starting point is the closed shell ground state of YbCH + 3 molecule to which one electron is added within the correlation procedure. In this work, we added the electron to the lowest energy σ orbital, thus using a minimal model space. The consideration of higher energy orbitals through a larger model space was not shown to have a significant effect on the value of W d in YbF and YbOH [37,38].
We apply equations (6) and (7)  both the W d and W s factors, using λ de = 10 −8 a.u. and λ ks = 10 −7 a.u. for BaCH 3 and λ de = λ ks = 10 −6 a.u. for YbCH 3 . These field strengths were selected guarantying numerical stability as it is shown in Appendix VI B. In all cases, the convergence criteria of the coupled cluster energy as well as the Hartree-Fock energy was fixed in the 1 × 10 −11 − 5 × 10 −11 a.u range.
In the following, we present the recommended W d and W s values and their uncertainties and discuss the scheme we use to determine the latter, focusing separately on the various parameters that determine the quality of the calculations. Initially, we focus our study on W d , and subsequently, we include W s in our discussion.

Electron correlation
In this work, we investigate the effect of various computational parameters within the relativistic coupled cluster approach on the calculated W d and W s factors.
a. Correlation space Table II and Figure 2 present the CCSD results obtained correlating a different number of electrons; in all these calculations, the virtual cut-off was adjusted symmetrically to the number of electrons correlated (that is the positive energy cut-off was taken to be of the same absolute size as the negative energy cut-off). We observe that correlating only the outer-corevalence electrons (17 electrons in BaCH 3 ; Ba: 5s; 5p; 6s, C: 2s; 2p, H: 1s, and 29 electrons in YbCH 3 ; Yb: 5p; 4f; 6s, C: 2s; 2p, H: 1s) causes relative errors of ∼9% compared to the value obtained when all the electrons are included in the correlation treatment. Previously, we found that correlating only outer-core-valence electrons in the isoelectronic BaF also led to an error of ∼10% compared to the all-electon calculation [37].
Therefore, for the recommended values, all the electrons were correlated. Note that the strong dependence on the description of the core region is unusual for most other atomic and molecular properties, such as molecular geometries or spectra.
Correlation of all the electrons in the coupled cluster approach requires simultaneous inclusion of a proportionally large number of virtual orbitals. Notice that the symmetric virtual cut-off is ∼ 1400 a.u. in BaCH 3 and ∼ 2300 a.u. in YbCH 3 . Figure 3 and Table III  the calculated W d of BaCH 3 at the CCSD and CCSD(T) levels using the 2z, 3z, and 4z quality basis sets (Table IV). In all the cases, the inclusion of triple excitations slightly reduces the calculated W d . The contribution of the triple excitations increases with the cardinality of the basis set, as expected, but even at the 4z basis set cardinality reaches only 1.3%. The magnitude of this contribution is close to what was observed in previous calculations on BaOH (1.4%) and BaF (1.5%) and suggests that the triple excitations can be neglected when the computational resources are an issue. We report the final W d and W s at the CCSD(T) level and use the triples contribution to estimate the uncertainty due to neglecting higher-order excitations. The use of the FSCC approach allows us to obtain accurate values of W d and W s in a open shell system like YbCH 3 . However, as perturbative triple excitations are not yet implemented in the FSCC module of the DIRAC program, we can not use this method to perform a reliable calculation of the contribution of the triple and higher-order excitations. Therefore, we do use the difference between the single reference CCSD and CCSD(T) results in BaCH 3 to estimate the relative uncertainty due to the neglect of the triple excitations in YbCH 3 molecule.

Basis sets
Next to the electron correlation, the quality of the basis set plays a crucial role in accurate theoretical description of molecular properties. In the following subsections we analyse the effect of the cardinality and the special features of the basis sets on the calculated W d and W s factors. a. Complete basis set limit Table V presents the dependence of the calculated W d constants on the cardinality of the basis set. Overall, W d changes monotonically with increase of basis set cardinality. In BaCH 3 it decreases by 6.5% when going from 2z to 4z cardinality basis set, while In YbCH 3 it increases by 4.1%. We extrapolate our results to the complete basis set (CBS) limit, using the usual three-point Dunning-Feller e −αn scheme (n = 2, 3, 4) [65,66] for extrapolating the DHF energies and the two-point Helgaker et al. n −3 scheme (n = 3, 4) [67] for extrapolating the correlation energies. We also tested the Martin (n + 1 2 ) −4 scheme [68] and the recent scheme of Lesiuk and Jeziorski (based on the application of the Riemann zeta function) [69] for extrapolating the correlation energies. The latter three extrapolation schemes (labeled H, M and L in Table V, respectively) give evenly spread and very similar CBS limits. We use the central Helgaker CBS limit for our final results and the respective 95% confidence interval (1.96σ) based on the spread of the three schemes as our extrapolation uncertainty estimate. The same methodology was used in our previous studies [70].
b. Core correlating and diffuse functions The results shown in Table VI showcase the effects of including outer-and inner-core correlating (cv3z, aev3z) and diffuse functions (aug-v3z) on the calculated W d factors. The difference between the cv3z/ae3z basis set and the v3z basis set lies in the addition of tight functions with high angular momentum, namely, 2/5 f and 1/2 g functions for Ba atom and 1/1 d function for C. On the other hand, the difference between the aug-v3z and the v3z lies in the even-tempered addition of diffuse functions for each angular momentum in all the atoms. For BaCH 3 the two augmentation schemes cause changes that are similar in magnitude but opposite in sign, leading to a negligible net effect. The W d in YbCH 3 shows a negligible depen- dence on both types of augmentation. In this case, the difference between the dyall.v3z and the dyall.cv3z basis sets is small since the two sets are identical for Yb and the only differences are in the number of the corecorrelating functions on carbon. In the dyall.ae3z basis set for Yb, only one extra g function is added to the dyall.cv3z basis. Due to the small (and opposite) effects, we proceed with the non-augmented dyall.v3z basis sets for both molecules and we included the effect of the augmentation and including correlating functions in our uncertainty estimation.

C. Recommended values and uncertainty estimation
The extensive computational study carried out in the previous section allows us to determine the most suitable method for obtaining the recommended values of the W d and the W s constants of the two molecules.
For BaCH 3 , we provide the final value of W d calculated at the CCSD(T) level of theory using the dyall.v3z basis set, correlating all the electrons and including virtual orbitals up to 2000 a.u. giving the base value of 3.33 × 10 24 hHz ecm . To this, we add the CBS extrapolation correction of −0.11 × 10 24 hHz ecm evaluated correlating 17 electrons with a virtual space cut-off of 30 a.u. We estimated the uncertainty due to the finite basis set (cardinality, core correlating and diffuse basis functions) and due to the neglect of higher-order excitations based on the study in the previous sections. To determine the size of each source of error, we use the difference in W d obtained with the final method and W d obtained with a lower approximation (for a given computational parameter) as is schematically shown in Table VII. In addition, to estimate the uncertainty stemming from the molecular geometry optimisation, we calculated the W d factor using the experimental geometry range (Ba-C bond distance = 2.557 -2.570Å [61]) and found that the maximum variation with respect to the value obtained using the optimised geometry is 1.08%. Finally, the Dirac-Coulomb Hamiltonian used in this work assumes an instantaneous (non-relativistic) electron-electron interaction. The Breit interaction is the first-order correction, and as it contains two-electron operators [71] it is not possible to separate this effect from the electron correlation effects. While aware of this limitation, we evaluate effect of the dominant Gaunt term of the Breit interaction [72] at the Dirac-Hartree-Fock (DHF) level of theory resulting in an decrease in W d of 1.8%. We thus set an estimate of the uncertainty on our recommended values due to the missing Breit interaction and higher order effects to 1.8%. All the individual contributions to the uncertainty are summarised in Table VII and we estimate the total uncertainty (assuming the effects are independent) at 3.44%. The recommended value for BaCH 3 is thus W d = 3.224 ± 0.111 × 10 24 hHz ecm . For YbCH 3 , we provide the final value of W d calculated at the FSCC level of theory, correlating all the electrons and including virtual orbitals up to 2000 a.u, as we did in BaCH 3 . In YbCH 3 however, it was computationally unfeasible to correlate all the electrons when using the dyall.v3z basis. Nevertheless, since the effect of the correlation space is bigger than the effect of the size of the basis set, we employed the dyall.v2z basis set to obtain the base value of 13.25 × 10 24 hHz ecm and to this we added the CBS correction of +0.56 × 10 24 hHz ecm , evaluated correlating 29 electrons with a virtual space cut-off of 10 a.u. We used an analogous scheme to that employed for BaCH 3 to estimate the uncertainty of this result. However, since the experimental geometry is yet not known for YbCH 3 , we use the same relative uncertainty due to the geometry optimisation as we derived for BaCH 3 . This assumption is justified by the fact that we used the same methodology to optimise the geometry of the two molecules. Furthermore, to estimate the effect of the excitation rank, we include in our uncertainty estimation twice the relative uncertainty obtained for BaCH 3 to account for neglecting both the triples and the higher-order excitations. The recommended value for YbCH 3 is thus W d = 13.799 ± 0.349 × 10 24 hHz ecm . Figure 4 presents the relative uncertainties for the different sources discussed above and calculated as described in Table VII. The highest contributions to the total uncertainty in BaCH 3 is due to the basis set incompleteness and the current limitation to describing the Breit interaction correctly. Notice that the addition of diffuse and extra correlating functions have an opposite effect, meaning that we over-estimate somewhat this source of uncertainty. In the case of YbCH 3 , the employed basis set is highly converged. Consequently, the total uncertainty in W d is smaller for YbCH 3 than for BaCH 3 with the leading source of uncertainty coming from the electron correlation description.
We have also performed calculations of the W s parameters using the computational scheme employed for the recommended W d values. The different computational parameters included in the uncertainty estimation in this work have effects of the same magnitude on the calculated W d and W s in BaF [37], with the Gaunt term being the only exception. In BaF, the Gaunt term, and therefore, the total uncertainty, is smaller in W s than in W d . Consequently, in this work we assume the relative uncertainty of the W s values is on the same order of magnitude as the uncertainty found for the corresponding W d factors. All the recommended values are summarised in Table VIII   similar molecules.

D. Comparison and molecular bond analysis
The final calculated enhancement factors in BaCH 3 and YbCH 3 are of the same order of magnitude as in their corresponding isoelectronic linear molecules BaOH, BaF, YbOH, and YbF, and the non-linear polyatomic molecules BaOCH 3 and YbOCH 3 , respectively. However, the enhancement factors in BaCH 3 and especially in YbCH 3 are larger than in the other molecules.
To investigate the origin of this difference, we calculated the W d and W s parameters for all the systems using the same approach and used the Quantum Theory of Atoms In Molecules (QTAIM) [75] and the Natural Bond Orbital (NBO) [76] analyses to study the relevant bonds. Table IX presents W d and W s obtained at the CCSD(T) and FSCC(0,1) level for the Ba-and Yb-containing molecules, respectively. To conserve computational effort, the dyall.v2z basis set was used and a cutoff was set at -2 to +10 a.u. for Ba-and -1 to +10 a.u. for Yb-containing systems. For the QTAIM and NBO analysis, we employ the SR-ZORA Hamiltonian with the PBE functional [77] and the QZ4P basis set in the Amsterdam density functional (ADF) package [78].
In the QTAIM analysis, intermolecular interactions can be characterized according to the topology of the electron density at the bond critical points (BCPs) denoted by r c . Specifically, the Espinoza and coworkers criteria [79] classify the bonding interactions according to the ratio of the potential V (r c ) and kinetic G(r c ) energy density at the BCPs, |V (r c )|/G(r c ). A ratios < 1 correspond to long-range, ionic or hydrogen bonds, values between 1 -2 to intermediate bonds (with both ionic and covalent character) and ratios > 2 to covalent bonds. The ratio |V (r c )|/G(r c ) in Table IX suggests that the M-X (M = Ba, Yb; X = C, O, F) bonds can not be characterized as purely ionic bonds, and especially, the M-C  It means that in BaCH 3 and YbCH 3 molecules, the unpaired electron that experiences the P, T-odd interaction is more localized on the heavy atom, leading to a higher enhancement factor than for the other Ba-, Yb-containing molecules. We tested this conjencture by including YbH in our analysis, which follows the expected trend -the less electronegative hydrogen leads to more covalent bonding and an increase in the W d and W s values. A similar observation was made in previous works [80,81].

IV. EXPERIMENTAL CONSIDERATIONS
BaCH 3 and YbCH 3 are both prime candidates for creation of intense beams via buffer gas cooling [82]. The rotational spectrum of the ground electronic state of BaCH 3 has been studied via millimeter/submillimeter absorption [83] and was created via the reaction of barium metal with Sn(CH 3 ) 4 . Other alkaline-earth monomethyl molecules have been created via reaction of metals with other monomethyl species [84], including chloromethane (CH 3 Cl) which could be used to react with ablated Yb or Ba metal in a cryogenic buffer gas cell. While YbCH 3 has not been studied spectroscopically, the chemical similarity of Yb with the alkaline-earth metals means that these production methods would likely work as well, as has been demonstrated with a number of other Yb-containing molecules [34,85].
The P, T-violation measurement with these species would be performed in the K = 1 rotational state of the ground electronic and vibrational state, in which K−doubling splits the rotational states into doublets of opposite parity [27,34,86]. The state with K = 1 corresponds to the molecule rotating about its C 3v symmetry axis with one quantum of angular momentum. Alkalineearth metals bound to −CH 3 and −OCH 3 ligands have this state at roughly 160 GHz above the ground state [86][87][88], giving them lifetimes estimated to be minutes or longer. As this energy corresponds to ≈ 8 K, there will be appreciable population in a buffer gas source at a few Kelvin. The doubling in small, open-shell molecules with a metal-centered electron and a CH 3 group is dominated by anisotropic dipole-dipole interactions between the H nuclear spins and unpaired electron spin [34,89]. For CaOCH 3 this has been measured to be around 300 kHz [89], but due to the closer distance between the metal and H 3 group in MCH 3 molecules this value is likely a few times larger, around 1 MHz. Note that the overall level structure, including the states which would be used for a spin precession measurement, are similar to other symmetric top molecules [86]. Also similar is that the EDM sensitivity would saturate to ≈ 50% of the maximum values at the small fields required to mix the parity doublets [86,90].

V. CONCLUSIONS
In this work, we report the enhancement factors W d and W s needed for the interpretation of possible eEDM measurements on the promising molecules BaCH 3 and YbCH 3 . We carried out a systematic study to devise a computationally feasible scheme that provides accurate predictions along with realistic uncertainties.
We showed that the scalar-relativistic level of theory in combination with the coupled cluster method and the ANO-RCC basis sets is a feasible and accurate approach for the optimisation of small polyatomic molecules. Correlating the valence and core-valence (n − 1) and (n − 2) electrons gives essentially the same results as correlating all electrons. The geometry obtained with this approach is in a very good agreement with the available experi-mental reports (for BaCH 3 ).
The values of W d and W s calculated in this work for both BaCH 3 and YbCH 3 are slightly larger than in other Ba-and Yb-containing molecules. The relation between the increased covalent character of the heavy atom bond and the size of the W d and W s factors provides important insight for the search for promising candidates for precision experiments.
Accurate calculations of properties needed for interpretation of precision measurements in polyatomic molecules is a challenging task and it is necessary to methodically evaluate the effects of the most relevant computational parameters. In this work, we provide such analysis for the BaCH 3 and YbCH 3 molecules. We show that calculating the W d and W s values requires an accurate description of the electron correlation. Specifically, in the coupled cluster approach, the number of electrons included in the correlation description plays the most significant role. The size of the basis set has a comparatively smaller but nonnegligible effect. The systematic evaluation of the effect of various computational parameters also allowed us to estimate the uncertainty in the values presented in this work.

Contracted vs. uncontracted basis set
In electronic structure calculations, contracted basis functions constructed from Gaussians perform for many properties with similar accuracy but with dramatically reduced computational costs compared to uncontracted basis sets. In the geometry optimisation of BaCH 3 and YbCH 3 (Table X), we observe that the scalar-relativistic approximation introduces a smaller error compared to the effect of the basis set contraction. The difference between the results obtained using contracted ANO-RCC.VDZ and uncontracted dyall.v2z basis sets at the SR level is 1.7%. However, this difference is reduced to 0.8% when comparing the larger SR-ANO-RCC.VTZ and X2C-dyall.v3z methods. This is expected, since with the increasing cardinality, the basis sets become more saturated and thus the negative effect of contraction should decrease. We thus proceed with the computationally less expensive ANO-RCC basis sets. Based on the results in Table X, when selecting basis sets for geometry optimisations considering their cost-benefit values, it is advisable to choose a higher cardinality contracted basis set over an uncontracted but lower cardinality one.

Electron correlation
Table XI contains comparison between the results obtained within the CCSD and the CCSD(T) approaches. The calculations were carried out within the singlereference framework and 37 electrons were correlated. We observe that inclusion of perturbative triple excitations has a minor effect on the obtained geometry (0.2% and 0.3% at the DZ and TZ basis set cardinality, respectively) and thus conclude that optimisation on the CCSD level is of sufficient quality when computational resources are of importance.
Furthermore, as expected, we find that the molecular geometry is a property mainly dependent on the valence region of the molecule.  (6) and (7), the W d and W s factors can be calculated as the first derivatives of the energy with respect to the corresponding perturbation. To calculate the derivative using numerical differentiation, it is necessary to determine the field strength at which the total energy depends linearly on the perturbation. Figure 5 shows the dependence of the total energy in BaCH 3 on the λ de perturbation. Linear behavior is found at smaller fields of the order λ de = 10 −8 a.u. Similarly, for the W s factor, the total energy is found to be linear at λ ks = 10 −7 a.u. In YbCH 3 , λ de = λ ks = 10 −6 a.u. correspond to the linear regime. To support these small fields, the convergence criterion of the coupled cluster energy as well as the Hartree-Fock energy was fixed at 10 −11 a.u.