$P,T$-odd effects in YbCu, YbAg, and YbAu

The current best description of elementary particles and their interactions is given by the Standard Model (SM) of particle physics.¹ This model is capable of describing almost all experimental observations (excluding gravitational interactions) and of accurately predicting a wide range of diverse phenomena, which is why over time it has been consolidated as a well-tested physical theory. However, it does not address some important observed effects, such as the matter–antimatter asymmetry in the universe, the neutrino oscillations, and the existence and nature of dark matter and dark energy.^2,3 Over the past decades, many new theories and extensions of the SM have been proposed to explain these phenomena.^3,4 Testing and restricting these theories is important for advancing our understanding of the fundamental laws of physics.

A promising way to test some of these theories is to search for effects due to the simultaneous non-conservation of spatial $(P)$ and time-reversal $(T)$ parities in atoms and molecules, such as those arising from the electric dipole moments (EDMs) of electrons.⁵ Interactions between these EDMs and electromagnetic fields violate both temporal and spatial invariances.⁶ Within a $CPT$-invariant theory (⁠$C$ refers to charge-conjugation symmetry), if $T$ symmetry is violated, then the combined $CP$ symmetry must also be non-conserved so that clearly $P,T$ violation implies $CP$ non-conservation.² 

The sources of $CP$ violation described by the SM lead to the prediction of a free-electron EDM d_e of ∼5.8 × 10⁻⁴⁰ e cm.⁷ However, more sources of $CP$ violation beyond those predicted by the SM are needed to explain, for example, the observed matter–antimatter asymmetry. The additional sources of $CP$ violation would in turn lead to an increase in the magnitude of the electron EDM (eEDM), even bringing it into the reach of present-day precision experiments.² Experimental searches for these phenomena are currently being carried out in atoms and molecules, taking advantage of the enhancement of the atomic and molecular EDMs. These arise from $P,T$-violating interactions, mainly those taking place between the eEDMs and the large internal atomic or molecular electric fields, and also other $CP$-odd nucleon–electron and nucleon–nucleon interactions.^2,8–10 Currently, the lowest upper limit of the eEDM is set at |d_e| < 2.1 × 10⁻²⁹ e cm. This upper limit was reported after measurements conducted by the National Institute of Science and Technology (NIST) on the HfF⁺ molecular ion¹¹ and combining these results with those obtained by the ACME Collaboration in their analysis of the ThO molecule.¹² It is important to stress that if $CP$ violation is assumed to arise exclusively from d_e in the NIST experiment (i.e., if the scalar–pseudoscalar nucleon–electron couplings are neglected), then the upper limit would be |d_e| < 4.1 × 10⁻³⁰ e cm. This lowest upper limit has already put considerable constraints on some of the theories beyond the SM.^13,14

In the ground state of a molecule with zero nuclear spins and a single unpaired electron, there are two main contributions to the energy arising from $P,T$-violating interactions, that is, the interactions between the EDMs of the electrons and the electromagnetic fields and from the $P,T$-odd scalar–pseudoscalar nucleon–electron (S-PS-ne) neutral-current interactions.^9,15 Since the effects of the eEDM and the S-PS-ne interactions are enhanced by the molecular electronic structure, the two corresponding molecular enhancement factors W_d (related to d_e) and W_s (which enhances the S-PS-ne interactions) are of particular interest.

The choice of molecule for the measurements has a significant impact on the sensitivity of the experiment due to, among others, the system-dependent enhancement of the $P,T$-odd effects. Since in open-shell molecules that contain only one heavy element, this enhancement scales roughly as the cube of the atomic number of the heavier nucleus to which the unpaired electron is strongly linked to Refs. 10, 16, and 17, some heavy-element-containing molecules have an advantage over other systems. Furthermore, practical experimental considerations play a crucial role in selecting a candidate for experiments. For example, the use of ultra-cold molecules increases the interaction times and, hence, the experimental precision.¹⁸ Therefore, the laser-coolability of the selected molecule provides a clear advantage. Various molecular properties relevant for precision measurements (e.g., laser coolability and sensitivity to the measured phenomena, but also many others) can be determined theoretically before experimental investigations in support of such experiments. In particular, the enhancement factors cannot be measured and must be provided based on accurate electronic structure calculations.

In this work, we investigate the sensitivity of the YbCu, YbAg, and YbAu molecules to $P,T$-violating phenomena. These systems are of particular interest since they contain two metal atoms. Theoretically, this means that these atoms can be laser-cooled separately and then associated into a molecule afterward, eliminating the need to laser-cool the molecule as a whole.^19–21 So far, laser-cooling of Yb and Ag atoms has been demonstrated.^22,23 Laser-cooling of Cu and Au has not been demonstrated yet, but cooling schemes have been proposed.²⁴ Furthermore, these polar molecules have large molecular-frame electric dipole moments, due to the large electronegativity of coinage-metal atoms, allowing for their easy polarization. In particular, YbAg is considered a promising candidate for a next-generation clock-transition eEDM measurement.²¹ 

No experiments have been performed so far on YbCu, YbAg, or YbAu, but several experimental groups pursue the ultra-cold formation of other Ag-containing molecules.^25–28 However, high-accuracy calculations of the potential energy curves, molecular-frame electric dipole moments, electric quadrupole moments, and static electric dipole polarizabilities of the present systems have recently been performed.²⁹ Here, we employ the four-component (4c) relativistic coupled-cluster (CC) approach to calculate the enhancement factors of the $P,T$-violating interactions between the eEDMs and the electric fields in the systems, W_d, and of the $P,T$-odd S-PS-ne interactions, W_s. The enhancement factors are determined for the ground states, $X2Σ12+$⁠, of the three molecules. We also carry out an extensive computational study to assign uncertainties to the calculated values.

In open-shell molecules containing at least one non-zero nuclear spin I, some internal nuclear interactions lead to nuclear spin-dependent molecular $P,T$-violating effects, such as magnetic interactions between the electrons and the nuclear magnetic quadrupole moment (NMQM), which appears for nuclear spin I ≥ 1.^30,31 These effects are, however, outside the scope of this work.

Section II covers the main theoretical aspects of this work, detailing how both enhancement factors can be obtained from effective Hamiltonians. Section III contains a description of the methods employed to calculate these factors, as well as the scheme used for geometry optimization. Next, Sec. IV presents the obtained enhancement factors and their dependence on effects such as the choice of the nuclear charge density model, the method for the treatment of electron correlation, the choice of the basis set, and the internuclear distances. Finally, Sec. V contains a concise summary of our findings.

To compute the enhancement factors, the effects of the eEDM and the S-PS-ne interactions are taken as perturbations on the 4c relativistic Dirac–Coulomb (DC) Hamiltonian. Since d_e and k_s,K are small quantities, the effects arising from both interactions are minute. Therefore, first-order perturbation treatment will already yield highly accurate results.

The effects on hydrogenic atoms arising from a $P,T$-odd interaction between the permanent EDM of a single electron, parallel to its spin, and an electromagnetic field have been studied for the first time by Salpeter in 1958.³⁵ In his seminal work, he introduced a $P,T$-odd perturbation term corresponding to a permanent eEDM into the one-electron Dirac equation in a Lorentz-covariant formulation. This term is analogous to the so-called Pauli moment interaction term (representing a QED interaction of the lowest order, for non-relativistic energies, between the electromagnetic field and the Pauli anomalous electric and magnetic dipole moments of the electron) but pre-multiplied by the pseudoscalar Dirac operator γ⁵ (see pp. 47–51 of Ref. 36).

In atoms or molecules with at least one electron whose spin is unpaired, the effects arising from these $P,T$-violating interactions produce atomic or molecular (permanent) $P,T$-violating electric dipole moments that may be significantly larger than that of a free electron. In the past, it has been shown that the $P,T$-violating EDM of these many-electron systems is mainly influenced by the electrostatic interactions between the eEDM of the unpaired electron and the internal electric fields.³⁷ Therefore, when calculating the molecular enhancement parameter W_d [see Eq. (1)], it is possible to ignore effects such as the interactions between the eEDM and the magnetic fields and also the electron–electron Breit interactions.

When neglecting both, the interactions of the eEDMs with internal and external magnetic fields and the Breit interactions, the mean value of the Salpeter Hamiltonian can be equated to the expectation value of two different operators. One of these two effective Hamiltonians is the sum of only one-electron operators, while the other one also contains two-body operators. These two-electron contributions, however, are considerably smaller than the one-electron ones in that particular effective Hamiltonian, and as such, they can usually be safely neglected.^8,37

Unlike the effective Hamiltonian $Ĥeff-1eEDM$ of Eq. (11), which explicitly depends on the nuclear charge density distributions ρ_K(r), the Hamiltonian $Ĥeff-2eEDM$ given in Eq. (13) does not have a functional dependence on effects due to finite nuclear sizes. However, when perturbation theory is applied to calculate the molecular enhancement factor W_d using S2, these effects are implicitly accounted for through their explicit inclusion in the unperturbed Hamiltonian of the system, $Ĥ(0)$⁠.

Furthermore, S2 only contains one-electron operators within an approximation in which both the magnetic interactions and the electron–electron Breit interactions are neglected, reducing the computational complexity. However, a drawback of this effective Hamiltonian is in the fact that the non-relativistic limit of its mean value is not zero, while that of the Salpeter Hamiltonian vanishes.³⁹ Furthermore, S2 does not allow for the analysis of separate nuclear contributions to the calculated W_d parameters. Meanwhile, within S1, the non-relativistic limit of W_d is zero. While S1 includes two-electron contributions [see Eq. (8)], these are much smaller than the corresponding one-electron counterparts^17,37 and can usually be neglected. Therefore, S1 allows us to study the effective contributions to W_d arising from each nucleus of the system, keeping the correct non-relativistic behavior, and using only one-electron operators.

The calculations reported in this work were performed using the DIRAC 19.0 program package.^41,42 The bond lengths were optimized using the exact two-component (X2C) Hamiltonian,⁴³ the single-reference coupled-cluster method with single and double excitations (CCSD), and the s-aug-v4z basis set.^44–46 The active space energy cutoffs for the virtual (unoccupied) and occupied orbitals were set to ±20 E_h, ±10 E_h, and ±10 E_h for YbCu, YbAg, and YbAu, respectively. A smaller active space was used for the heavier molecules since computing their bond lengths is computationally more intensive. The optimized internuclear distances for YbCu, YbAg, and YbAu are 2.7543, 2.8589, and 2.6524 Å, respectively. The bond length increases from YbCu to YbAg and then decreases from YbAg to YbAu due to relativistic effects. A non-relativistic treatment would show longer bonds for heavier systems.⁴⁷ The enhancement factors were computed for these calculated equilibrium internuclear distances.

The molecular enhancement factors W_d and W_s were calculated using the four-component DC Hamiltonian and the multi-reference Fock-space coupled-cluster method with single and double excitations (FSCC).⁴⁸ A multi-reference method is required for these calculations due to the challenging character of the ytterbium-containing molecules.^49,50 This method has been employed previously to study the enhancement factors of the $P,T$-violating interactions in YbOH⁵¹ and YbCH₃.⁵² In addition, the FSCC implementation within the EXP-T program^53,54 was used to investigate the effects of including triple excitations.

The uncontracted Dyall’s valence-only basis sets of double-ζ (v2z), triple-ζ (v3z), and quadruple-ζ (v4z) quality were employed.^55,56 In addition, the core–valence basis sets (cvXz, X = 2, 3, 4) were used to examine effects arising from the correlation of core electrons.⁵⁷ These particular basis sets add tight functions with large exponents to the vXz basis sets. Finally, the augmented basis sets (s-aug-vXz) were employed to investigate how accurately the outer regions of the systems were described. The s-aug-vXz basis sets add a single diffuse function to each symmetry block.

For calculations of W_s, identical nuclear models (either PN or GN) were used in both the $Ĥ(0)$ and $ĤS-PS-ne$ operators [see Eqs. (2) and (19)]. In addition, when W_d was calculated using S1, the PN model was employed in both $Ĥ(0)$ and $Ĥeff-1(PN)eEDM$ [see Eqs. (2) and (12)]. The GN models were not yet implemented in DIRAC within $Ĥeff-1eEDM$⁠. Finally, computations of W_d within S2 were performed using either PN or GN models in the unperturbed Hamiltonian $Ĥ(0)$⁠.

A field strength λ = 10⁻⁶ was set when studying both parameters W_d and W_s in YbCu and YbAg, whereas λ = 10⁻⁷ was used for the calculations involving YbAu.

To obtain accurate enhancement factors accompanied by well-defined uncertainties, multiple computational aspects will be addressed. For all the calculations discussed in this section, the 4c DC Hamiltonian was employed. First, the baseline results will be presented in Sec. IV A. Then, the effect of selecting two different nuclear charge density distribution models is examined for both factors in Sec. IV B. Furthermore, the two schemes employed for computing W_d will be compared in Sec. IV C. Next, the influence of the basis set will be determined in Sec. IV D. Thereafter, different computational approaches will be compared in Sec. IV E. Finally, the effect of the geometry of the system on W_d and W_s will be discussed in Sec. IV F.

This extensive investigation allows us to set uncertainties on the final values. The justification for the final results and their uncertainties will be given in Sec. IV G.

All the reference values for W_d and W_s were computed on the v3z/FSCCSD level, using the equilibrium bond lengths obtained in this work and Gaussian-type nuclear charge distribution models. The reference values for W_d have been obtained using S2 [see Eq. (13)]. For YbCu and YbAg, a virtual space cutoff of 500 E_h was employed, and two and four electrons were frozen, respectively. For YbAu, the virtual space cutoff was set to 40 E_h, and 56 electrons were frozen. The selection of these correlation parameters is justified in Table S1 of the supplementary material.

The calculated enhancement factors for the three systems are given in Table I. We can observe that both parameters are very similar for YbCu and YbAg, while much lower absolute values are obtained for YbAu. This finding will be elucidated in Sec. IV C. Furthermore, multiple corrections to these baseline values will be determined. The final obtained values will be given in Sec. IV G.

The effects arising from the use of two different nuclear models in the calculations of the molecular enhancement factors were investigated. For all three molecules, these effects were analyzed in both W_d and W_s. In this study, S2 has been employed to calculate W_d. All calculations were performed on the FSCCSD/v2z level of theory, correlating all electrons and using a virtual space cutoff of 3000 E_h for all three systems.

In Fig. 1, the effects of using different nuclear charge density distribution models on W_d (S2) are shown (see also Table S3 of the supplementary material). For YbCu and YbAg, the results using GN are $∼1.1%$ smaller than those obtained employing PN, while for YbAu, the GN result is $∼9.2%$ larger than the PN one. In particular, the nuclear size effects on YbCu and YbAg agree with those on YbF previously reported by Gaul and Berger in Ref. 17. Employing S1, they modeled the internal electric field of Eq. (11) using PN and GN models (using in both cases GN to model the unperturbed Hamiltonian) and found that these two results for YbF differ by $∼0.9%$⁠.

Meanwhile, in Fig. 2, the nuclear size effects on W_s are shown. The contributions from each of the nuclei are shown, as well as the total values of W_s, represented by hatched blocks. The contributions associated with the nuclei of the coinage metals become increasingly important as their atomic number increases. Moreover, as expected, these contributions are of opposite sign to those of ytterbium because the gradients of the nuclear densities have opposite directions in the region between the two nuclei. For YbAu, the two contributions almost cancel each other out, leading to a very small total value of W_s (see Table S4 of the supplementary material).

In all three cases, the contribution of the ytterbium nucleus is reduced by $∼13%$ when going from PN to GN, while for copper, silver, and gold, the reduction is ∼0.16%, $∼3.5%$⁠, and $∼20%$⁠, respectively. Thus, the effect of using a finite nucleus model becomes more significant for the heavier elements, as can be expected. The total absolute W_s is also lower for the calculations performed using the GN model.

It is well known that when the electric field produced only by the nuclei is taken into account in S1, the two schemes described in Sec. II A should yield similar results for W_d.³⁷ While the use of S2 is computationally less demanding, since it requires just a single calculation per system instead of the two that are required for each diatomic molecule, S1 allows us to examine the effective contributions arising from each nucleus. To study those individual contributions, we present (to the best of our knowledge) the first four-component results of W_d using the approximate effective Hamiltonian of Eq. (12), corresponding to the use of S1.

The molecular enhancement factors W_d were computed using both schemes (S1 and S2) on the FSCCSD/v2z level of theory, using the PN model. Symmetric cutoffs of 500 E_h, 500 E_h, and 40 E_h were employed for YbCu, YbAg, and YbAu (freezing 2, 4, and 56 electrons), respectively. It is worth mentioning that the differences between the values of W_d computed according to S1 and S2 (using GN for the unperturbed Hamiltonian and the internal electric field of the perturbation in S1) have been previously studied using two-component methods in a set of almost 30 diatomic molecules.¹⁷ 

In Fig. 3, we show the differences in W_d when computed with the two schemes for the three systems (see also Table S5 of the supplementary material). The values of W_d for YbCu and YbAg calculated using S1 are about 1.5% higher than those obtained using S2; this difference is 5% for YbAu. These differences may arise from contributions due to two-electron interactions, which are not taken into account in our implementation of S1 but are implicitly included in S2, but also from our use of finite basis sets, which could lead to different representation qualities of the different involved operators. These contributions are on the order of 1% for molecules containing one predominant heavy element (YbCu and YbAg), in agreement with previous predictions.^8,17,37,60

Using S1, it can be seen that the contributions to W_d from the Yb nucleus remain almost constant for all three systems. Moreover, the contributions from the coinage metals all have opposite signs to those coming from Yb, as is the case for the W_s factors. This is expected since the electric fields due to the two nuclei in the internuclear region have opposite directions, and according to Eqs. (11) and (12), this generates opposite contributions to this enhancement parameter. The decreasing total W_d factor from YbCu to YbAu is due to the increasing contribution from the second nucleus (opposite to that of the first).

To observe how the size of the basis set influences the enhancement factors, W_d (S2) and W_s were computed with double-ζ, triple-ζ, and quadruple-ζ quality basis sets. These calculations were done using GN models, FSCCSD method, and the (occupied and virtual) active space cutoffs were set to ±20 E_h, ±10 E_h, and ±10 E_h, freezing 38, 64, and 82 electrons of YbCu, YbAg, and YbAu, respectively.

The plot in Fig. 4 shows the effect of increasing the basis set cardinality on W_d (see Table S6 of the supplementary material). Deviations of at most 3.5% are observed when this enhancement factor is calculated for YbCu and YbAg using v2z and v4z basis sets, whereas these differences reduce to less than 0.5% when v3z and v4z results are compared. While no apparent convergence can be observed for the total W_s values of YbAu, this convergence can be seen by looking at the individual contributions in Fig. 5. A similar trend is expected for W_d.

Apart from adding contributing functions to all orbitals, it is also possible to add only tight or diffuse functions. Tight functions should increase the accuracy of the description of the core region of the system, while diffuse functions should improve the accuracy of the description of the valence region of the system. For these calculations, we used the cvXz basis sets, which contain higher angular momentum tight functions, and the s-aug-vXz basis sets, which augment a diffuse function per symmetry to the vXz basis set. These computations were done using the same computational settings as the ones above.

The effects arising from increased accuracy in describing the core and valence regions of the systems are small for both enhancement factors. Adding tight functions has a negligible effect on all three systems for W_d (with the largest change of 0.4%) and W_s (with the largest change of 0.3%). Meanwhile, adding diffuse functions reduces the parameter values slightly (with the largest changes of 2.9% and 1.9% for W_d and W_s, respectively). Furthermore, the enhancement factors of the YbAg system are more affected than those of YbCu. For YbAu, the addition of diffuse functions increases the calculated W_d and decreases the W_s. Still, the relative effects remain within 3% of the total value, as can be seen in Table S6 and Fig. S4 of the supplementary material.

For the final calculations, the v3z basis set was used for baseline results, and then, corrections due to basis set quality and a better description of the core and valence regions were added. The incompleteness of the basis sets was also taken into account in the treatment of uncertainties.

The method used so far throughout this work for calculations of the enhancement factors was 4c Dirac–Coulomb FSCC. Other approaches to calculate W_d (S2) and W_s have been analyzed, and the results are compared in Table II and Figs. 6 and 7. These methods are DHF, Møller–Plesset up to second order (MP2), FSCCSD, and FSCCSDT. The calculations were performed employing the dyall.v2z basis set and GN models, and the active space cutoffs for the post-DHF computations were set to ±2 E_h, freezing 66, 84, and 116 electrons of YbCu, YbAg, and YbAu, respectively. This small active space and basis set were used because the inclusion of triple excitations (using the EXP-T program) is computationally expensive. The DHF results differ the most from the FSCC values, while the MP2 values are close to the FSCC results for YbCu and YbAg. Meanwhile, going from MP2 to FSCC reduces the value of W_d for YbAu by a factor of two and reverses the sign of W_s.

Including the triple excitations in the FSCC calculations has only a minor effect on the enhancement factors of YbCu and YbAg, but increases the W_d and decreases the W_s of YbAu significantly, in line with the enhanced sensitivity of these systems to the other computational parameters. It should be pointed out that due to the small size of the basis set, the contribution of the triple excitations should be considered as an order-of-magnitude estimate, rather than an accurate prediction. The differences between the FSCCSD and FSCCSDT methods will be used to estimate the uncertainty due to the neglect of the higher excitations in Sec. IV G.

The equilibrium bond lengths used in all the calculations of the enhancement factors reported so far were obtained in Sec. III, as no experimental bond lengths are available for the systems considered in this work.

To investigate the influence of the bond lengths on the enhancement factors, these parameters were computed at the calculated equilibrium internuclear distances reported in Sec. III (using S2 for W_d), as well as at 0.05 and 0.1 Å larger and smaller internuclear distances. These computations were performed on the FSCCSD/v2z/GN level of theory, and the (occupied and virtual) active space cutoffs were set to ±100 E_h, ±100 E_h, and ±95 E_h, freezing 12, 20, and 28 electrons of YbCu, YbAg, and YbAu, respectively. The results are given in Table III. The absolute values of both enhancement factors of YbCu and YbAg decrease from the values at equilibrium bond lengths by less than 1.8% for the larger internuclear distance and increase by 2.2% at most for the smaller bond length. For YbAu, the deviations are significantly larger, up to 64%.

Given the strong (and almost linear) dependence of the enhancement factors on the internuclear distances (see Table III), the bond lengths of YbCu, YbAg, and YbAu were re-optimized so that some effects that were not accounted for previously (especially perturbative triple excitations and a larger active space) can be considered. Since there are no experimental internuclear distances for our systems, we first analyzed six molecules for which there are available measurements with which our optimizations can be compared: YbX (X = F, Cl, Br)⁶¹ and YAu (Y = Be, Mg, Ca).^62–64 These optimizations were performed using the Molpro package,^65,66 using a spin-free X2C Hamiltonian, ANO-RCC basis sets (with an effective cardinality larger than quintuple-ζ),^67–70 and treating electron correlation effects by using the single-reference coupled-cluster method with a full treatment of single and double excitations and perturbative treatment of triple excitations [CCSD(T)], correlating all virtual orbitals and valence electronic orbitals (with a cutoff of −55 a.u.). The use of a scalar-relativistic framework is justified here as these molecules are studied in the $X2Σ12+$ states. Our tests showed that the differences due to the neglect of spin–orbit effects are below 0.01 Å, i.e., smaller than the effects of basis set size and electron correlation, which we were trying to address.

In all cases, these optimized bond lengths differ by less than 0.03 Å with respect to the available experimental values, with a mean absolute error of 0.012 Å (as can be seen in Table S11 of the supplementary material). Therefore, we can assume that the new internuclear distances obtained for YbCu, YbAg, and YbAu would have a similar error range. The same method was used to re-optimize the molecules of interest arriving at the final bond lengths of 2.7229, 2.8329, and 2.6975 Å for YbCu, YbAg, and YbAu, respectively. By fitting the enhancement factors and internuclear distances given in Table III with a linear function, we have corrected the values of W_d and W_s considering the new optimized distances. The results of these corrections are shown in Sec. IV G.

The anharmonicity observed in the potential energy curves for the ground electronic states of YbCu, YbAg, and YbAu implies that their effective equilibrium bond lengths are slightly shifted compared to the minimum of the potential energy curves. These slight differences in bond lengths result in small changes in the enhancement factors.

To compute the vibrational corrections to W_d and W_s, calculations of these enhancement factors and the potential energies were performed for different bond lengths. The results for W_s, split up into contributions arising from each nucleus, can be found in Fig. 8 (and also in Table S8 of the supplementary material). These values are given as a function of the difference between the internuclear distances and the equilibrium bond length of the corresponding molecule. The resulting W_s,Yb/X values are shown as a percentage compared to the results found at the equilibrium bond lengths given in Sec. III. It can be seen that the coinage metal contributions are significantly more sensitive to the bond-length effects compared to the Yb contributions. This is also in agreement with the sensitivity of the respective atoms to the description of the electronic structure, as can be seen in Fig. 7 (and also in Table S7 of the supplementary material).

The vibrational corrections were computed using the VIBCAL module available in DIRAC-19.0, using a fourth-order polynomial for the energy fitting, and employing enhancement factors calculated at the equilibrium bond length given in Sec. III, 0.01, 0.05, and 0.1 Å larger and smaller internuclear distances. Scheme 2 was used to compute W_d. The calculations were performed on the FSCCSD/v2z/GN level of theory, and the active space cutoffs were set to ±100 E_h, ±100 E_h, and ±95 E_h for YbCu, YbAg, and YbAu, respectively. The vibrational corrections are listed in Table IV for both enhancement factors, and they alter the baseline values of the enhancement factor (obtained at R_eq) by $<1%$⁠.

The baseline results given in Sec. IV A can now be corrected to obtain the final values of W_d and W_s, which can also be considered as recommended values based on this work. The contributions due to using a larger basis set, correlating all electrons, increasing the virtual space cutoffs, including triple excitations in the FSCC method, and taking into account bond length and vibrational corrections were added to the reference baseline values. Table V and Table S9 of the supplementary material present the values of these different contributions and the final values.

The basis set corrections were calculated as differences between FSCCSD calculations with active space cutoffs of ±20 E_h for YbCu and ±10 E_h for YbAg and YbAu, employing the dyall.v4z and dyall.v3z basis sets (see Table S6 of the supplementary material).

The effects arising from correlating all electrons were taken as the differences between freezing 2, 4, and 56 electrons for YbCu, YbAg, and YbAu respectively, and correlating all the electrons and simultaneously increasing the virtual cutoff to 6000 E_h at the FSCCSD/v2z level of theory. A detailed analysis of the effects of the active correlation space size is also available in Sec. I of the supplementary material.

The corrections due to the inclusion of higher excitations were taken as the differences between FSCCSDT and FSCCSD calculations, employing v2z basis sets and energy cutoffs of ±2 E_h. Finally, the effects arising from the corrected bond lengths were taken as described in Sec. IV F 1, and the vibrational effects were extracted from Table IV.

To provide conservative and reliable uncertainties for these enhancement factors, a treatment similar to that given in Refs. 52, 71, and 72 was employed. The uncertainties obtained from the discussed considerations are given in Table VI and Table S10 of the supplementary material and analyzed individually in Subsections IV G 1–IV G 5.

The uncertainties arising from the basis sets used in the calculations have three different sources: (i) the general quality of the basis sets (related to their cardinality), (ii) the number of extra tight functions, and (iii) the number of diffuse functions employed.

The overall quality of the basis sets affects the uncertainty of the final results because they cannot be considered fully converged at the v4z level. The baseline values are corrected by adding the differences between the v4z and v3z results, as described in Sec. IV D, and also adding the differences between the s-aug-v3z and v3z results and between the cv3z and v3z values.

For both parameters, the uncertainties in basis set quality are seen to mainly arise from the differences between the v4z and v3z results. These differences are taken as the uncertainties related to the basis set cardinality.

The uncertainties due to the possible insufficient amount of tight functions were taken into account as the differences between the cv3z and v3z results. Finally, the differences between the s-aug-v3z and v3z values were taken as the uncertainties in the treatment of diffuse functions.

Electron correlation is affected by the chosen (occupied and virtual) active space cutoffs, along with the rank of excitations taken into account. To account for the frozen orbitals in the baseline calculations, the results provided in Table S1 and Figs. S1 and S2 of the supplementary material are used. The differences in enhancement factors between correlating all electrons and freezing 2, 4, and 56 electrons for YbCu, YbAg, and YbAu, respectively, were used to correct the baseline values, as shown in Table V; we thus do not assign uncertainties related to not correlating all the electrons.

The corrections for the limited virtual space cutoff that is used in the baseline calculations are taken as the differences between the obtained values at virtual space cutoffs of 6000 E_h and 500 E_h for YbCu and YbAg and between 6000 E_h and 40 E_h for YbAu. The v2z basis set was used in these calculations, and the results are shown in Table V and also in Table S2 and Fig. S3 of the supplementary material.

To account for the effects of higher-lying virtual orbitals, half of the differences between the results obtained using virtual space cutoffs of 6000 E_h and 3000 E_h were taken as additional uncertainties (see Table VI).

The baseline values do not take into account triple excitations. The differences between the results obtained using FSCCSDT and FSCCSD methods, as described in Sec. IV E, were used to correct this omission. Since these results may not be fully converged yet, the effect of including triple excitations is used to estimate the uncertainty arising from not including quadruple and higher excitations, by multiplying by one-half the difference between the FSCCSDT and FSCCSD results.

Inaccuracies in the calculated values of internuclear distances contribute to the enhancement factors’ uncertainties. Based on our benchmark geometry calculations discussed in Sec. IV F 1, we conservatively set the uncertainty in the bond length to 0.05 Å and calculate the enhancement factors for internuclear distances ±0.05 Å around R_eq (see Table III). The resulting differences were taken as the uncertainties associated with the molecular geometry.

Since the interaction constants k_s,K in Eq. (19) are specific to each nucleus, calculating the total W_s value as a sum of the atomic contribution is only approximate and introduces an associated error. We analyze this error in detail in an upcoming publication.⁷³ Here, we provide the resulting uncertainties of 0.059 h kHz, 0.308 h kHz, and 1.313 h kHz for YbCu, YbAg, and YbAu, respectively. We note that the specific errors differ for each isotopologue. Therefore, here, we use the isotopic average weighted by the natural abundances of all constituting elements.

To compute the total uncertainties for the reference values of the enhancement factors of YbCu, YbAg, and YbAu, the Euclidean norm of the individual uncertainties is taken. These total uncertainties are obtained on the assumption that the different contributions are largely independent since they concern high-order effects.

The relative uncertainties for W_s are similar to those found for W_d. For these systems, the enhancement factors are dominated by the uncertainty on the bond lengths, followed by the uncertainties due to missing higher CC excitations and, to a smaller degree, by the basis set incompleteness. The calculated values of W_s and W_d of YbCu and YbAg have uncertainties of about a single percent. A similar computational approach yielded uncertainties of 2%–7% for the enhancement factors of BaF,⁷¹ BaCH₃, YbCH₃,⁵² and LaO, LaS, and LuO.⁷² 

The results for YbAu have significantly higher relative uncertainties than those of the other two systems. Both enhancement factors of this system are relatively small due to a cancellation of large similar-sized contributions from the two constituent atoms, rendering them very sensitive to computational settings and leading to large relative uncertainty.

The calculated enhancement factors are compared to those found for other ytterbium-containing molecules in Table VII. In addition, some of the systems currently used in experiments aiming to restrict the lowest upper limit on the eEDM were added.^11,12,74,75 The W_d and W_s factors of YbCu and YbAg are of similar magnitude to those found for YbOH, YbCH₃, YbOCH₃, and YbF. Both enhancement factors of YbAu are significantly smaller than those found for any other ytterbium-containing systems.

Since both the interactions of the eEDMs with electromagnetic fields and the S-PS-ne neutral-current interactions may contribute to an eventual experimental detection of $P,T$-violating effects in molecules, these two types of interactions should be decoupled from each other. This can be done by performing measurements on systems with different enhancement factor ratios.¹⁷ For the systems studied in this work, and some other molecules currently and previously under investigation, these ratios can be calculated from the values of enhancement factors reported in Table VII.

Interpretation of experiments that search for $P,T$-violating effects arising from both the interactions between eEDMs and electromagnetic fields and the S-PS-ne neutral-current interactions in molecules requires knowledge of the W_d and W_s enhancement factors, which can only be obtained through molecular electronic structure computations. In this work, these parameters were computed for the YbCu, YbAg, and YbAu systems, selected due to their possible experimental advantages. The enhancement factors W_d were calculated using two different schemes. This study confirms that, as indicated by previous predictions and studies,^8,17,37,60 the two-body contributions to W_d represent only about 1% of its value in the case of systems with a predominance of one heavy element (YbCu and YbAg) but reach up to 5% for YbAu, formed from two heavy elements. In addition, a thorough uncertainty analysis was performed to assign a conservative error to the obtained results.

The final values were calculated using the FSCC method and the 4c DC Hamiltonian in conjunction with relativistic basis sets. The main contributing sources of uncertainty are due to the limited basis set size and the neglect of CC excitations beyond triples.

The obtained enhancement factors of YbCu and YbAg are of very similar size to other Yb-containing compounds investigated in the literature. In the case of YbAu, the cancellation of the contributions arising from the two nuclei in the system leads to vanishingly small total W_d and W_s values. For YbCu and YbAg, the results are also of similar size as for other systems currently investigated experimentally to search for signs of $P,T$-violating effects. Compared to YbF and YbOH, the alternate method of producing and cooling these systems provides an alternative route for future experiments setting a new lowest upper limit on the eEDM.

In the supplementary material, we analyze the influence of the active space (in particular, of the energy of occupied and virtual correlated orbitals) on the calculation of Wd and Ws. We also present a set of tables where we report the values of these parameters using different nuclear models, different schemes corresponding to the use of the two effective Hamiltonians described in this work (for the case of Wd), and a few different basis sets. In addition, tables are provided showing the dependence of the molecular enhancement factors on the use of different methods to treat electron correlation and vibrational effects and also the contributions to Wd and Ws associated with each nucleus of the studied molecular systems.

We acknowledge the University of Groningen’s Center for Information Technology and the Dutch National Supercomputer for their support and for providing access to the Hábrók and Snellius high-performance computing clusters. This work made use of the Dutch national e-infrastructure with the support of the SURF Cooperative using Grant Nos. EINF-5787, EINF-8014, and EINF-8532. We thank the anonymous referees who have made high-quality comments and suggestions that have significantly improved this article. I.A.A. thanks R. Berger and K. Gaul for inspiring discussions and acknowledges the partial support from FONCYT through Grant Nos. PICT-2021-I-A-0933 and PICT-2020-SerieA-0052 and CONICET through Grant No. PIBAA-2022-0125CO. The work of I.A.A., A.B., and S.H. was supported by the project Probing Particle Physics with Polyatomic molecules with Project No. OCENW.M.21.098 of the research program M2, which is financed by the Dutch Research Council (NWO). The work of A.B. was also supported by the project High Sector Fock space coupled cluster method: benchmark accuracy across the periodic table with Project No. Vi.Vidi.192.088 of the research program Vidi, which is financed by the Dutch Research Council (NWO). The work of L.F.P. and S.H. was supported by the project Searching for missing antimatter with trapped molecules with Project No. VI.C.212.016 of the research program Vici, which is financed by the Dutch Research Council (NWO). L.F.P. acknowledges the support from the Slovak Research and Development Agency Project Nos. APVV-20-0098 and APVV-20-0127.

The authors have no conflicts to disclose.

Johan David Polet: Data curation (supporting); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Visualization (lead); Writing – original draft (lead); Writing – review & editing (equal). Yuly Chamorro: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Supervision (equal); Validation (equal); Writing – original draft (supporting); Writing – review & editing (equal). Lukáš F. Pašteka: Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Resources (equal); Software (equal); Validation (equal); Writing – review & editing (equal). Steven Hoekstra: Funding acquisition (equal); Investigation (supporting); Project administration (supporting); Resources (equal); Validation (supporting); Writing – review & editing (equal). Michał Tomza: Formal analysis (equal); Investigation (supporting); Methodology (equal); Writing – review & editing (equal). Anastasia Borschevsky: Formal analysis (equal); Funding acquisition (lead); Investigation (equal); Project administration (lead); Resources (equal); Supervision (equal); Writing – original draft (supporting); Writing – review & editing (equal). I. Agustín Aucar: Conceptualization (lead); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (lead); Methodology (equal); Resources (equal); Software (lead); Supervision (lead); Validation (equal); Visualization (supporting); Writing – original draft (equal); Writing – review & editing (lead).

The data that support the findings of this study are available from the corresponding author upon reasonable request.