The Role of Relativistic Many-Body Theory in Electron Electric Dipole Moment Searches Using Cold Molecules

: In this review article, we survey some of our results pertaining to the search for the electric dipole moment of the electron (eEDM), using heavy polar molecules. In particular, we focus on the relativistic coupled cluster method (RCCM) and its applications to eEDM searches in YbF, HgX (X = F, Cl, Br, and I), BaF, HgA (A = Li, Na, and K), and YbOH. Our results are presented in a systematic manner, by ﬁrst introducing the eEDM and its measurement using molecules, the importance of relativistic many-body theory, and ﬁnally our results, followed by future prospects.


Introduction
The electric dipole moment of the electron (eEDM), which is a parity (P) and time (T) reversal violating property [1,2], is of immense interest in probing fundamental physics, especially beyond the Standard Model (BSM) theories. The Standard Model (SM) predicts a small value for eEDM (∼ 10 −38 e-cm), as compared to many BSM theories, which predict the eEDM to be several orders of magnitude larger [3] (see Figure 1). Therefore, by comparing experimental upper bounds on the eEDM (the eEDM has not been detected yet, in over six decades of ever-improving experiments) with those predicted by several BSM theories, one can stringently constrain the latter. The general interest in this field is compounded by the fact that high precision eEDM search experiments are non-accelerator approaches that complement accelerator searches, but can explore much higher energy scales, as high as PeV [4]. The eEDM can also throw light on the baryon asymmetry in the Universe [5], which attempts to address why our Universe is matter-dominated, while one may expect an early Universe to be created with equal amounts of matter and anti-matter. This connection is associated with the fact that both the not-yet-measured eEDM and the well known value of baryon asymmetry coefficient are CP violating (T violation implies CP violation, due to the CPT theorem [6]). Further details can be found in Ref. [7].
The current best limits for the eEDM are from heavy polar molecules (all at 90 percent confidence intervals): ThO (d e < 1.1 × 10 −29 e-cm; d e = 4.3 ± 3.1 stat ± 2.6 syst × 10 −30 e-cm; the subscripts 'stat' and 'sys' refer to the uncertainties due to statistics and systematics, respectively) [8], followed by HfF + (d e < 1.3 × 10 −28 e-cm; d e = 0.9 ± 7.7 stat ± 1.7 syst × 10 −29 e-cm) [9], and YbF (d e < 1.06 × 10 −27 e-cm; d e = −2.4 ± 5.7 stat ± 1.7 syst × 10 −28 e-cm) [10,11]. Two experiments that employ BaF to probe the eEDM are currently underway [12][13][14]. An eEDM experiment that commenced recently (called the polyEDM experiment) using YbOH is expected to improve the sensitivity further in the future [15].  It turns out that there are no contributions at the two-and three-loop levels either; a sample three-loop diagram has been given above, and although an individual three-loop diagram itself may not vanish, the non-zero diagrams eventually cancel out [16], (c) The first non-zero contribution to the SM eEDM occurs at the four-loop level [17], and (d) supersymmetry (SUSY) contribution to the eEDM. The phases at the two vertices are different, and therefore, one can obtain a non-zero contribution to the eEDM at the one-loop level itself. This is also why the eEDMs in most beyond the SM (BSM) are much larger than their SM counterpart [3].
This article is structured as follows: since the introduction section covers the eEDM, we shall discuss about eEDM searches in molecules in the next section, followed by a brief discussion on relativistic many-body theory with focus on the relativistic coupled cluster method (RCCM), after which we shall present some of our results for YbF, HgX (X = F, Cl, Br, and I), BaF, HgA (A = Li, Na, and K), and YbOH. We shall use atomic units in this article, unless explicitly stated otherwise.

The Effective Electric Field
A high precision experiment that directly uses an electron beam, for example, would be extremely difficult, as the particle would accelerate away in the presence of an electric field. Instead, one uses atoms or molecules as proxies, and measures the shift in their energy level due to the eEDM. We shall only consider heavy polar molecules in this article, as molecules are preferred over atoms due to their higher experimental sensitivity. The shift in the energy of a molecule (∆E) that is prepared in some state, due to the eEDM (d e ), is given by where E eff is called the effective electric field, since it has the same dimensions as an electric field. The expression for E eff can be obtained by starting with the eEDM Hamiltonian [18] The above expression is basically the well-known -d · E expression for the interaction potential for a dipole, d, in an electric field, E, except that the dipole here is the eEDM that lies along the spin axis [19] (d e Σ i ; Σ i is the 4 × 4 Pauli matrix), while the field in this case is an internal electric field, E intl i . The internal electric field, E intl i , experienced by an electron, is due to the nuclei and all the other electrons in the molecule. The β Dirac matrix denotes that we work in a relativistic regime.
A relativistic treatment is necessary, since otherwise, the effective electric field is zero, due to the Schiff's theorem [20]. The summation is over the number of electrons, N e . Because the eEDM is an extremely small quantity, the shift in energy due to it, to first order in perturbation, can be expressed as In the first equation given above, |Ψ is the unperturbed wave function of a molecule. We note that although the eEDM operator is parity-odd, the resulting expression for ∆E with respect to the unperturbed wave function is non-zero, unlike an atomic case. This is because one expands the molecular orbitals (MOs) contained in a molecular wave function as a linear combination of 'atomic orbitals' (AOs; strictly speaking, basis functions), often called the Linear Combination of Atomic Orbitals-Molecular Orbital (LCAO-MO) approach. This leads to mixing of opposite parity orbitals, and therefore a non-zero contribution. From Equations (1) and (3), one obtains As seen earlier, the internal electric field, E intl i , experienced by the ith electron, includes contributions from each of the nuclei as well as from the other electrons. Since the latter leads to a term containing a two-body operator, directly calculating it is computationally very expensive. Hence, one employs an effective eEDM Hamiltonian,Ĥ eff eEDM , whose expectation value is same as that of the originalĤ eEDM , when Ψ is the exact eigen function of the unperturbed Hamiltonian. Further details can be found in Refs. [21,22]. Hence, the effective electric field is given by In the above expression, γ 5 is the product of Dirac matrices, while p i is the momentum operator.

eEDM: A Combination of Experiment and Theory
As seen earlier in Equation (1), an eEDM experiment measures a shift in energy ∆E, due to the eEDM, of a molecule prepared in a particular state. However, we need to revisit Equation (1) and slightly modify it, in order to accommodate for the dependence of the sensitivity to an external electric field that is applied to polarize the molecules. Therefore, one obtains We see from the above expression that η, which is called the polarizing factor, carries with it the information pertaining to the applied electric field (E extnl ), the molecular permanent electric dipole moment (also known as the PDM, denoted by D), and the energy difference between rotational states of the molecule, ∆E rot . Often, for 2 Σ systems, one defines a polarizing field, E pol = 2B/D (where B is the rotational constant of the molecule, and is related to ∆E rot ), which sets the scale for the required external electric field in order to polarize the molecules. It is desirable to be able to polarize a molecule with as low an external electric field as possible, or as low an E pol as possible. Therefore, a molecule with a very low PDM can actually hamper the prospects of an otherwise good eEDM experiment, and hence, polar molecules are the preferred choices. The shift in energy is measured in experiment, the polarizing factor can either be calculated or found from experiment, while E eff can only be calculated using a relativistic many-body theory. Therefore, an eEDM experiment is one of those branches in physics where experiment and theory go hand in hand. Also, it is worth noting at this point that the eEDM itself is a particle physics property that answers deep questions such as the matter-antimatter asymmetry in the universe (which is on a significantly different length scale), while the property can be measured by a combination of high-precision molecular experiments and relativistic many-body theory (involving quantum chemistry).

Choosing a Molecular Candidate
The statistical sensitivity of an eEDM experiment is given by N refers to the number of molecules used per coherence time. T refers to the total interrogation time, while τ is the coherence time of the state of interest to an eEDM experiment. The importance of E eff can be immediately understood from the above expression; it occurs outside the square root, and therefore plays a significant role in the sensitivity of an eEDM experiment.
The choice of a good eEDM candidate molecule usually relies upon the above mentioned factors (besides the effects of systematics and the prospects that a molecule offers for rejection of certain systematic effects). Not all of them can be satisfied in any one system; however, one looks for a suitable combination of these parameters to propose a new candidate which can potentially set a better bound on the eEDM. The effective electric field in a reasonably good molecular eEDM candidate is extremely large, and is usually ∼GV/cm (as compared to atoms, whose 'E eff ' for a ∼10 kV/cm field is ∼MV/cm). This is one of the major reasons why molecules are preferred over atoms.

Relativistic Many-Body Theory
We need to evaluate the expression given in Equation (5), in order to obtain E eff . The molecular wave function in the expression has to be obtained from a relativistic many-body theory. We employ the relativistic coupled cluster method (RCCM) for our accurate calculations. Solving the RCCM equations rely upon the relativistic Hartree-Fock (HF) or Dirac-Fock (DF) reference state, and therefore, we shall briefly explain the DF equations, and then motivate RCCM by drawing parallels with the more intuitive many-body perturbation theory (MBPT).

The Dirac-Fock Method
The HF equations and its relativistic counterpart, the DF equations, can be obtained by applying the variational principle, which guarantees that the energy thus obtained would be an upper bound to the actual ground state energy (for example, Ref. [23]). The equations arê There are N such equations for an N-electron system, each corresponding to the ith spin-orbital, |ϕ i (1) (the many-body wave function itself is built from these spin-orbitals, or single particle wave functions, as a Slater Determinant, |Φ 0 = Det{ϕ 1 ϕ 2 · · · }). The operator,t(1) contains in it the kinetic and the electron-nucleus interaction terms, andv(1, 2) is the familiar electron-electron interaction term. The summation j is over the number of electrons in the system. The second and third terms are called the direct and the exchange terms respectively. Written in this form, the Hd or DF equations reveal the physical meaning of the direct term. By rearranging the direct term, one can readily see that it leads to a mean field interpretation, that is, each electron experiences an average potential due to the remaining (N − 1) electrons. The exchange term has no classical counterpart. Also, in this form, it is transparent that the equation has to be solved in a self-consistent way, starting with an initial guess for the spin-orbitals.
In the DF method, one employs a relativistic Hamiltonian, that is, the Dirac-Coulomb Hamiltonian (DCB),Ĥ 0 , given byĤ The summation, i (and j), is over the electrons, while A is over the nuclei. α and β are the Dirac matrices, while p i is the momentum operator. Z A is the atomic number of the Ath nucleus. r i is the position vector from the origin to the site of the ith electron, and r A is that from the origin to the site of the Ath nucleus. In practice, the third term would correspond to a finite nucleus, and we have given the point nucleus case above for simplicity. We work in the Born-Oppenheimer approximation, where the nuclei are 'clamped', and also in the no-pair approximation [24], where one projects out the negative energy solutions from the electronic Hamiltonian. Also, the single particle wave function, for a relativistic case, is built from four-component spinors, given by where v and w are referred to as the large and small components respectively, and in practice, are related to one another by the kinetic balance condition. The condition guarantees that the expression for kinetic energy is recovered in the non-relativistic limit, and also aids in avoiding variational collapse [25,26].
We shall now briefly discuss the matrix formulation of the DF equations, by first rewriting Equation (8) (we drop the co-ordinate labels for simplicity) as followŝ Here, |ϕ i is a general spin-orbital. Letf ≡t +v DF . Then, Expanding |ϕ i as a linear combination of a suitable basis set, that is, |ϕ i = ∑ n C in |χ n , and projecting from the left with χ m |, one obtains.
The DF coefficients, C in , are what one solves for. S mn is the overlap integral. It is now straightforward to recast the above equations in a matrix form as After performing the similarity transform shown above, the resulting equation, that is, Equation (18), is solved iteratively to obtain the coefficients. Note that the choice of basis is one of the factors that determines the precision of a calculation. The other important factor that determines the precision is electron correlations, which we shall discuss next.

Electron Correlation
The HF method (and its relativistic counterpart, the DF method) is a mean field approach, which ignores the fact that every electron actually interacts with each of the remaining electrons. The physical effects beyond those embodied in the mean field approximation are subsumed under the umbrella of electron correlation.
If E is the energy that one obtains by solving the Schroedinger's equation exactly, and if E DF is the DF energy for the same system, then the correlation energy, E corr , is Some well known examples of many-body theories that take into account electron correlation effects are configuration interaction (CI), multi-configuration Hartree-Fock (MCHF), many-body perturbation theory (MBPT), coupled cluster Methods (CCM), and their relativistic extensions.

Many-Body Perturbation Theory: Non-Relativistic and Relativistic
To introduce MBPT, we define zeroth-order Hamiltonian (Ĥ 0 ) as a sum of Fock operatorŝ For this Hamiltonian, the Hartree-Fock determinant (Φ 0 ) and its k-particle k-hole excitation determinants (Φ k ) become eigen functions ofĤ 0 , and hencê In this picture, a hole refers to an occupied orbital, while a particle to an unoccupied/virtual one. Electron correlation effects are taken into account by accommodating for such virtual excitations to the HF wave function. Also, E k = E 0 − ∆E k ; where the second term contains the orbital energies associated with a k particle-k hole excitation. For example , for one-particle one-hole and two particle-two hole cases, we obtain respectively. In this notation, i, j, k, · · · are holes, and a, b, c, · · · are particles.
We now introduce a perturbation to the Hamiltonian, which is the difference between the actual 1 |r i −r j | and the HF potential. Then, the Schroedinger's equation corresponding to the actual wave function of the system would beĤ where |Φ (1) 0 contains one order ofĤ , |Φ (2) 0 contains two orders ofĤ , etc. This framework is common for both the relativistic (denoted by the subscript, 'rel') and non-relativistic cases (denoted by the subscript, 'NR'), with the orbitals being different (four-component), and the one-body part of the Hamiltonian also differing as follows withV Nuclear i being a shorthand for ∑ A Z A |r i −r A | for a point nucleus, or a suitable expression for a finite nucleus, depending on whether it is a Fermi or a Gaussian one.

Relativistic Coupled Cluster Method
The RCCM is a very accurate approach to treat electron correlation effects. To that end, the method has been referred to as the gold-standard of many-body theory of atoms and molecules. The RCCM wave function is given by T is called the cluster operator. For an N-electron system,T =T 1 +T 2 + · · · +T N . Here, T 1 = ∑ i,a t a iâ †î , where all possible one-particle one-hole excitations are accounted for by the summation, t a i is called the t-amplitude, which is like a weight factor for a one-particle one-hole excitation, andâ †î is a shorthand notation forâ † aâi , where a hole 'i' is annihilated and a particle 'a' is created. Similarly, T 2 = ∑ i =j,a =b t ab ijâ †b †ĵî , etc. One must also include all possible combinations of two independent one-particle one-hole excitations (given byT 2 1 ), a one-particle one-hole and a two particle-two hole occurring independently, but again all such possibilities (that is,T 1T2 ), and so on. By taking into account all such options and accounting for over-counting, we naturally obtain an exponential structure, that is, Equation (30) [27].

RCCM and MBPT
We shall now compare RCCM with MBPT. We begin by examining the full MBPT wave function, and then connecting it with RCCM.
Equations (32) and (33) involve rewriting the corrections to the wave function at a given order of perturbation, as a summation over single (that is, one-particle one-hole, and denoted by subscript S), double (D), · · · excitations. If we plug in Equations (32) and (33), and the rest of the corrections to the wave function to infinite orders, and rearrange, the resulting equation as We see that this is actually the coupled cluster wave function. This correspondence is allowed, as the wave functions obtained from full MBPT and CCM without any truncation have to be the same. Equation (34) shows that for a given level of particle-hole excitation, CCM contains in it those excitations to all orders of perturbation in the residual Coulomb interaction. RCCM is also both size extensive and size consistent. Further details can be found in Ref. [28].

The Energy and Amplitude Equations
The RCCM equations that are to be solved, in order to obtain the t amplitudes, can be derived as followsĤ and, Φ * 0 |e −TĤ eT|Φ 0 = 0.
In Equation (38), we have projected Equation (37) with the DF wave function and in Equation (39), |î †ĵ †bâ , and so on. One can also write e −TĤ eT as (ĤeT) c , according to the linked cluster theorem. We shall briefly outline the derivation of the linked cluster theorem here. Given an operatorÂ, e −TÂ eT can be expressed as a sum of nested commutators, by using the Baker-Campbell-Hausdorff expansion. By using suitable arguments from the second quantization formalism, one can show that each of the terms in the sum are fully contracted, or in the diagrammatic approach, fully connected (denoted by the subscript, 'c'), leading to the final expression e −TÂ eT = (ÂeT) c [29]. Finally, we note that the expression for energy is naturally terminating, as this shall have important consequences later in this article.

The CCSD Approximation
One often works with the coupled cluster singles and doubles approximation (CCSD, or its relativistic counterpart, RCCSD), whereT =T 1 +T 2 . In such a case, the energy and amplitude equations are After solving the amplitude equations and obtaining the wave function of the system, the next step would be to obtain the property of interest, for example, E eff . This can be done in two ways, by either solving an expectation value problem, or by adopting an energy derivative approach. We shall discuss both these approaches below.

The Expectation Value Approach
The expectation value of an operator,Â, in RCCM is given by The last equation was derived by Cizek [30,31]. The subscript 'N' means that the operator is normal ordered [32], while 'C' stands for fully connected, as seen earlier. This approach has the advantage of making the underlying physics more transparent, as one can break an expectation value expression down into its constituent terms, which in turn can be associated with several many-body effects, such as Bruckner pair correlation. However, the expectation value expression given above is an infinite, non-terminating series, and therefore requires that one manually truncate the expression, for obtaining a result. One such approximation is the linear expectation value approximation (LECC), which we have adopted. The expectation value in the CCSD approximation is then given by As we shall see later, this turns out to be a good approximation for computing E eff , for many systems with one unpaired electron.

The Energy Derivative Approach
The second approach relies on evaluating a property as an energy derivative. The derivative is taken with respect to a perturbation parameter, λ. The total Hamiltonian,Ĥ(λ), is given bŷ Then,Ĥ (λ)|Ψ(λ) = E(λ)|Ψ(λ) ; Also, From perturbation theory, Within the RCCM framework, this approach takes advantage of the fact that the expression for coupled cluster energy naturally terminates. A first order property like E eff or PDM could be evaluated as an energy derivative, either analytically (i.e., analytical derivative) or numerically (i.e., the finite field method). In our work, we have implemented the latter. We note that our discussions are based on computing first order properties like E eff and PDM. The approach can be extended in a straightforward way to higher order properties, such as polarizabilities and hyperpolarizabilities.
We shall now explain some details of our implementation. We use an unrelaxed version of finite field coupled cluster method (FFCC), that is, we solve the unperturbed DF equations, and add the perturbation at the CCSD level. We shall derive an expression for CCSD correlation energy (E corr ) with the perturbation lambda, as follows.Ĥ

(λ)|Ψ(λ) = E|Ψ(λ)
Using the fact that Φ 0 |Ĥ|Φ 0 = E DF , and operating with e −T on both sides on the left, we obtain Using the linked cluster theorem, our correlation energy is We note that the effects associated with orbital relaxation would be recovered at the perturbed CCSD level; for example, see Ref. [33]. The next step is to evaluate the energy derivative. Depending on the degree to which we wish to minimize the truncation error in a calculation, we can opt for a two-point central difference formula or beyond to evaluate the derivative. The two-point and the six-point formula are, respectively The property of interest is obtained by adding its DF contribution to the above formula, that is, The above equation can be obtained for a forward difference formula, for example, by starting with E(λ) and E, and taking the difference. If the perturbation is d e , andÂ isĤ eff eEDM /d e , then one obtains E eff , or with an electric field as the perturbation, withÂ being the dipole operator, one gets the PDM.
This approach offers the advantage of not having to truncate an infinite series in order to obtain a property. However, it is extremely time consuming, as one has to carefully choose the right value of λ, and a two-point central difference formula means that one has to evaluate the CCSD equations twice for one property.

Methodology
For all our results, we performed the DF calculations and obtained the property integrals using the UTChem package [34][35][36]. The same program was also used for atomic orbital to molecular orbital integral transformation. We then employed the Dirac08 program [37] to solve the CCSD equations and obtain the t amplitudes, and finally, we used our expectation value code for the LECC computations. For a finite field calculation, we introduce the perturbation at the CCSD level. Therefore, the one-electron integrals as well as the relevant property integrals with a chosen value of the perturbation parameter are fed to the Dirac08 program to obtain the energies. Further details on LECC approach can be found in Ref. [21], while Ref. [38] provides more information on the FFCC approach.

First Application of RCCM to YbF
YbF currently holds the third best limit on the eEDM, and improvement in its current sensitivity is expected in the near future. We performed our first accurate molecular RCCM calculations on YbF, and using the LECC approach, obtained a value of −23.1 GV/cm for E eff (see Table 1), with an estimated error of less than 10 percent [21]. A fairly large quadruple zeta (QZ) basis (primitive Gaussian basis sets) was used (Yb: (35s,30p,19d,13f,5g,3h,2i), F: (13s,10p,4d,3f)) [39][40][41]. We used a bond length of 2.0161 Angstroms [42]. Our result was an improvement over the previous ones, as they were based on the DF approximation [43], MBPT [44], effective core potential methods [45], and the CI method in the singles and doubles approximation [46]. The latest result on the molecule uses a complex generalized Hartree-Fock and a complex generalized Kohn-Sham scheme within the zeroth order regular approximation [47]. We had also systematically investigated the behaviour of E eff as well as the PDM with change in the number of virtuals and frozen orbitals. Additionally, we also presented the individual correlation contributions, from Equation (46), at the QZ level of basis. In conclusion, our accurate result leads to an increase in the value of the upper bound on the eEDM to 11.8 × 10 −28 e-cm, as compared to the earlier result of 10.5 × 10 −28 e-cm. Further details can be found in Ref. [21].

Mercury Halides
It was seen from Equation (7) that a large value of E eff would substantially improve the statistical sensitivity of an eEDM experiment. We had identified that mercury halides, HgX, possess E eff ∼ −100 GV/cm (see Table 2) [48], which is one and a half times larger than that of the current best candidate, ThO, and about five times that of YbF and HfF + . We estimate the error in our calculations to be within 5 percent. We employed the Dyall c2v basis sets for Hg [49], and cc-pVDZ (correlation consistent-polarized valence double zeta) basis for the halides [50], for our computations. We had chosen the following bond lengths (in Angstroms): HgF (2.00686) [51], HgCl (2.42), HgBr (2.62), and HgI (2.81) [52]. Our HgF calculation was an improvement over those by Meyer et al. [53,54], and Yu Yu Dmitriev et al. [55]. These systems also possess favourable experimental features, and ultracold HgX can be produced in large quantities, for example, by photoassociation of laser-cooled Hg with magnetically trapped halogen atoms [56]. Anticipating N ≥ 10 3 molecules/s, and τ ∼ 1 s with a slow beam or fountain of HgX molecules (for example, Ref. [57]), the eEDM sensitivity, δd e ≤ 3 × 10 −27 / √ N e-cm. We also found that E pol is less than 10 kV/cm for HgX, and thus one could apply small lab electric fields to polarize the molecules. We also point out that HgX molecules have sets of internal comagnetometer states, which can lead to suppression of systematic effects. In a recent work, Yang et al. explore laser cooling possibilities for HgF, and estimate the statistical sensitivities for an eEDM experiment to be ∼9 × 10 −31 e-cm (beam experiment) and ∼10 −32 e-cm (trap experiment) [58]. Further details on our calculations and estimates can be found in Ref. [48]. Table 2 also presents our results for HgX, using FFCC [38]. We employed both the two-and the six-point formula, and the results did match perfectly. We also chose the following values of λ: 10 −3 , 10 −4 , 10 −5 , and 10 −6 , and found a smooth convergence in the results. The results show that LECC is a good approximation, for calculating E eff . We took HgF as a representative case for obtaining our error estimates in our calculations, and broadly categorized the errors under three sources: 1. percentage fraction difference between the TZ (triple zeta) and DZ (double zeta) FFCC results was 3.02, 2. percentage fraction difference between DZ cases with and without diffuse functions was 2.4, and 3. percentage fraction difference between FFCCSD(T) and FFCCSD results, at DZ level of basis was 3.3. Table 2. Summary of our calculated LECC and finite field coupled cluster method (FFCC) values of E eff for HgX and BaF. We used double zeta (DZ) basis sets for HgX, while BaF's result is using a QZ basis, with diffuse functions added to Ba. All the units are in GV/cm.

BaF
As mentioned earlier, BaF is a promising candidate for an eEDM search experiment, and an experiment using the system is underway [12]. Also, a second experiment on the same molecule has commenced, which plans on trapping a large number of BaF molecules in a solid rare-gas matrix. This approach is expected to drastically improve the statistical sensitivity of an eEDM experiment [13,14]. We employed the FFCC method, and obtained −6.50 and −6.46 GV/cm, using LECC and FFCC respectively (see Table 2) [38]. We used QZ basis sets, from Dyall (and appended diffuse functions to it from Sapporo basis) [41,59] and BSE databases (cc-pV basis sets) [50], for Ba and F respectively. We used a bond length of 2.16 Angstroms [60,61]. Our calculations were an improvement over the earlier ones by Kozlov et al. [62], Nayak and Chaudhuri [63], and Meyer et al. [53,54].

Mercury Alkalis
We performed LECC calculations on E eff of mercury alkalis (HgA; A = Li, Na, and K), and found that their effective fields are comparable with that of YbF's (see Table 3) [64]. We used Dyall's cVTZ basis for all of the molecules [65,66]. The bond lengths (in Angstroms) were: HgLi: 2.92, HgNa: 3.52, and HgK: 3.90 [67]. In combination with a reasonably large E eff , we found based on a preliminary survey that HgA satisfy the criteria of very long τ and a large N. An eEDM experiment with a long τ can be envisioned with ultracold molecules trapped in optical lattices [68]. We conservatively estimate that N ∼ 10 4 ultracold HgA molecules can be produced in an optical lattice, based on what was demonstrated with isoelectronic YbLi molecules [69,70]. We also find that the values of E pol are (71,28,17) kV/cm for (HgLi, HgNa, HgK), which implies that it is feasible to polarize trapped, ultracold, HgA molecules. With an integration time of ∼10 7 seconds, a preliminary set of estimates for the sensitivity, δd e = (1.3; 2.5; 3.1) × 10 −30 e-cm for (HgLi, HgNa, HgK) respectively. Further details on our calculations and estimates can be found in Ref. [64].

YbOH
Polyatomic molecules are currently emerging as promising eEDM candidates. The advent of polyatomic systems was marked by the first fast sisyphus laser cooling of SrOH on the experimental front [71], followed by a proposal using RaOH for an eEDM experiment [72]. A proposal by Kozyryev and Hutzler, followed immediately by commencement of the experiment, brought YbOH to the forefront of eEDM searches [15]. In their work, the authors highlight the importance of YbOH; the possibility of laser cooling as well as possessing bending modes with closely-spaced parity doublets. The implications due to the former is that laser-cooled YbOH molecules in an optical lattice would tremendously improve the coherence time to ∼ 1s on the statistical sensitivity front, while drastically reducing systematics, such as motional magnetic fields. The latter implies internal co-magnetometer states, which eliminates the systematics associated with reversing electric fields. Lastly but very importantly, the spectroscopy of the system is already known. These factors lead the authors to expect a sensitivity that is four orders better than the current best candidate, ThO. However, the sensitivity also relies upon a reasonably large E eff . In our work on YbOH, we calculate E eff values not only for a linear geometry, but also for a bent one, as the state of interest is a bent mode. Table 4 shows the obtained results in a linear geometry, with Dyall's double (DZ), triple (TZ), and quadruple zeta (QZ) basis, with the Yb-O bond length being 2.0026 Angstroms, and O-H 0.922 Angstroms [73]. The result is similar to that of YbF. Also, there were two other works on linear geometry, at about the same time as ours on YbOH: Gaul and Berger (who used a zeroth order regular approximation) [47,74], and Denis et al. (who also employ a relativistic coupled cluster method) [75]. For an excited bent mode, specifically a (010) mode, we do not know beforehand the degree to which the molecule is bent, and therefore, we estimate the effective electric field at the DF level for various values of Θ, which is the angle formed by O-Yb, with respect to H-O, as Figure 2 shows. A DF treatment is sufficient for our purposes, as the dominant contribution to E eff comes from it. We assume that the correlation contributions would not vary significantly with change in angle. Contrary to expectation, our results show that the DF effective electric field could vary significantly in a low-lying bent mode, between about 13 and 18 GV/cm, and also shows an oscillatory behaviour, for the chosen points. For further details, refer to [76].

Future Prospects
We shall now look at two of the many-body methods that we plan on implementing in the near future, due to their advantages over the conventional coupled cluster approach.

The Normal Coupled Cluster Method
In this variant of the coupled cluster method, the bra and the ket are defined on different footings [77], asĤ where |Ψ = eT|Φ 0 ;T = ∑ I =0 Here, T I s refer to the t amplitudes, andĈ + I to the operators associated with particle-hole excitations.Ĉ − I are the set of creation and annihilation operators, which are associated with de-excitations.

Amplitude Equations by Variational Principle
We can derive the amplitude equations for the normal coupled cluster method (NCCM), starting from the variational principle. Consider an energy functional, given by According to the variational principle: Using ∑ I =0 T IĈ − I =ˆ T − 1, we obtain: In the second term on the left hand side,Ĉ + I and e −T can be switched, since they commute. Therefore,

The Expectation Value
Since the bra and the ket are treated differently, as seen earlier, the resulting NCCM expectation value expression naturally terminates. For any operator,Â, this can be seen as follows: For illustration, let us consider NCCM in the singles and doubles approximation, that is, NCCSD. One can show that Therefore, e −TÂ eT terminates depending on the rank ofÂ, and therefore NCCM expectation value avoids the issue of truncation.

The Hellmann-Feynman Theorem
An important attribute of NCCM that makes it attractive, besides being variational and having an expectation value that naturally terminates, is that it satisfies the Hellman-Feynman Theorem (HFT), that is, whereĤ(λ) is the Hamiltonian depending on a continuous parameter λ, |Ψ(λ) is an eigen state of the Hamiltonian, and E λ is the energy of the state |Ψ . For NCCM, the Hellmann-Feynman theorem takes the following form: This expression can be rigorously proved. ExpandingĤ(λ), Ψ(λ) and Ψ(λ) in the above equation, withĤ(λ) =Ĥ 0 + λÂ, and equating the zeroeth powers of λ, we obtain Here, Ψ and Ψ are the zeroth order or unperturbed parts of the wave functions for the bra and ket in lambda CC respectively. Therefore, the effective electric field could either be calculated as an expectation value, as shown earlier, or as an energy derivative, as shown above. For example, in a recent work, Sahoo and Das calculate the enhancement factors of the Hg atom using spherical symmetry to probe nuclear EDMs in the expectation value approach [78], while Sasmal et al. compute the effective electric field of RaF using the energy derivative approach [79]. We are planning on implementing the former and the extended coupled cluster method (ECCM) for our molecular calculations. In the ECCM, the same ket is used as in the NCCM, but theˆ T in the bra is replaced with eˆ T [77].
Starting from here, one can show that after some algebra and equating terms of the order lambda, the energy and amplitude equations can be obtained as Therefore, one solves the coupled cluster equations, hence obtaining the t amplitudes, followed by solving for the amplitude equations given above to obtainT (1) . Finally, from the energy equation, one can evaluate E 1 . The effective electric field could thus be obtained, as it is a first order property.
Having introduced the two methods, we shall now compare them with the conventional coupled cluster approach. As seen earlier, the NCCM expectation value naturally terminates, unlike the CCM expression. Therefore, one expects NCCM to be more accurate (for a given level of particle-hole excitation and choice of basis), as it can capture missing correlation effects in CCM that arise as a consequence of manually terminating an infinite series. NCCM also satisfies the Hellmann-Feynman theorem. The analytical derivative approach recasts a property as an energy, and therefore there is no truncation involved in this method either, hence making this approach more accurate than CCM. We also note that the analytical derivative approach fares better than FFCC, since the latter involves running several calculations to obtain a property and also could have numerical issues. However, both NCCM and analytical derivative are harder to implement, and are computationally more expensive, as both require us to solve for an additional set of amplitudes. The NCCM treats the bra and the ket differently, and we need to evaluate the unperturbed bra and ket amplitudes, while in the analytical derivative approach, we need to solve for unperturbed and perturbed ket amplitudes. In CCM, one solves only for the unperturbed ket amplitudes. Also, using these methods for polyatomic systems is expected to be more challenging, and the degree of complexity could vary depending on factors such as the state of the molecule, and the number of valence electrons in that system.

Conclusions and Outlook
Relativistic many-body theory is indispensable in the searches of electron EDM. It is necessary to determine the upper limit for eEDM and also for identifying promising candidates for eEDM searches. With new results expected for the electron EDM experiments, improvements in relativistic many-body calculations in molecules are desirable. Relativistic coupled cluster theory is well suited for the calculations associated with eEDM searches mentioned above. It would be necessary to develop different variants of the theory for this purpose: normal CC method, the analytical derivative approach, ECCM [77,81], and tensor network tailored CC theory for open-shell molecules [82].
Author Contributions: All of the authors have contributed equally in the execution of the reported work. V.S.P. prepared the manuscript, B.P.D. supervised this review article, and the other co-authors were involved in discussing the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.