Zeroth order regular approximation approach to electric dipole moment interactions of the electron

TL;DR

This paper introduces a quasi-relativistic two-component computational approach using ZORA for efficient calculation of electron electric dipole moment interactions, validated against more complex relativistic methods on heavy diatomic molecules.

Contribution

It presents a novel two-component method combining cGHF and cGKS within ZORA for accurate and efficient eEDM interaction calculations, reducing computational complexity.

Findings

Method accurately reproduces relativistic four-component results.

Purely relativistic effects are captured by the upper component treatment.

Validated on molecules like RaF, BaF, YbF, and HgF.

Abstract

A quasi-relativistic two-component approach for an efficient calculation of $\mathcal{P,T}$-odd interactions caused by a permanent electric dipole moment of the electron (eEDM) is presented. The approach uses a (two-component) complex generalized Hartree-Fock (cGHF) and a complex generalized Kohn-Sham (cGKS) scheme within the zeroth order regular approximation (ZORA). In applications to select heavy-elemental polar diatomic molecular radicals, which are promising candidates for an eEDM experiment, the method is compared to relativistic four-component electron-correlation calculations and confirms values for the effective electrical field acting on the unpaired electron for RaF, BaF, YbF and HgF. The calculations show that purely relativistic effects, involving only the lower component of the Dirac bi-spinor, are well described by treating only the upper component explicitly.

Tables

Table 1. Table 1: Comparison of literature data of the 𝒫 , 𝒯 𝒫 𝒯 \mathcal{P,T} -odd eEDM interaction parameter W d subscript 𝑊 d W_{\text{d}} of the spin-rotational Hamiltonian of 137 BaF calculated with different four-component methods and with a quasi-relativistic GHF/GKS-ZORA approach.

Method	$\frac{W_{\text{d}}e\cdot \text{cm}}{{10}^{24}\text{Hz}\cdot h}$
Exp.+SE1   (Ref. Knight, 1971; Kozlov and Labzowsky, 1995)		$-3.5$
Exp.+SE1   (Ref. Ryzlewicz et al., 1982; Kozlov and Labzowsky, 1995)		$-4.1$
GRECP/SCF2a   (Ref. Kozlov et al., 1997)		$-2.30$
GRECP/SCF/EO2b   (Ref. Kozlov et al., 1997)		$-3.75$
GRECP/RASSCF2c   (Ref. Kozlov et al., 1997)		$-2.24$
GRECP/RASSCF/EO2d   (Ref. Kozlov et al., 1997)		$-3.64$
DF3a   (Ref. Nayak and Chaudhuri, 2006b)		$-2.93$
DF-RASCI3b   (Ref. Nayak and Chaudhuri, 2006b)		$-3.52$
NR-MRCI4   (Ref. Meyer, Bohn, and Deskevich, 2006)		$-2.5$
	$W_{\text{d},\text{I}}$		$W_{\text{d},\text{II}}$
GHF-ZORA (this work)	$-3.3$		$-3.3$
GKS-ZORA/B3LYP (this work)	$-2.9$		$-2.9$

Table 2. Table 2: Comparison of literature data of the 𝒫 , 𝒯 𝒫 𝒯 \mathcal{P,T} -odd eEDM interaction parameter W d subscript 𝑊 d W_{\text{d}} of the spin-rotational Hamiltonian of RaF calculated with different four-component methods and with a quasi-relativistic GHF/GKS-ZORA approach.

Method	$\frac{W_{\text{d}}e\cdot \text{cm}}{{10}^{24}\text{Hz}\cdot h}$
SODCI1   (Ref. Kudashov et al., 2014)		$-24.0$
FS-RCCSD+${\Delta }_{\text{basis}}$+${\Delta }_{\text{triples}}$2  (Ref. Kudashov et al., 2014)		$-25.6$
DF-CCSD3   (Ref. Sasmal et al., 2016b)		$-25.4$
	$W_{\text{d},\text{I}}$		$W_{\text{d},\text{II}}$
GHF-ZORA (this work)	$-28.0$		$-27.3$
GKS-ZORA/B3LYP (this work)	$-25.1$		$-24.4$

Table 3. Table 3: Comparison of literature data of the 𝒫 , 𝒯 𝒫 𝒯 \mathcal{P,T} -odd eEDM interaction parameter W d subscript 𝑊 d W_{\text{d}} of the spin-rotational Hamiltonian of YbF calculated with different methods and with a quasi-relativistic GHF/GKS-ZORA approach.

Method	$\frac{W_{\text{d}}e\cdot \text{cm}}{{10}^{24}\text{Hz}\cdot h}$
Exp.+SE+corr.1  (Ref. Van Zee et al., 1978; Kozlov and Ezhov, 1994; Kozlov et al., 1997)		$-12.6$
GRECP/SCF2a   (Ref. Titov, Mosyagin, and Ezhov, 1996)		$-9.1$
GRECP/RASSCF2b   (Ref. Titov, Mosyagin, and Ezhov, 1996)		$-9.1$
GRECP/RASSCF/EO2c   (Ref. Mosyagin, Kozlov, and Titov, 1998)		$-12.06$
GRECP/RASSCF/EO+4$f$2d    (Ref. Mosyagin, Kozlov, and Titov, 1998)		$-12.16$
RDHF+CP3    (Ref. Quiney, Skaane, and Grant, 1998)		$-6.0$
UDF4    (Ref. Parpia, 1998)		$-12.03$
DF5a    (Ref. Nayak and Chaudhuri, 2006a)		$-9.63$
DF-RASCI 5b    (Ref. Nayak and Chaudhuri, 2006a)		$-10.88$
DF-MBPT(2)5c    (Ref. Nayak and Chaudhuri, 2007)		$-10.43$
DF-RASCI* 5d    (Ref. Nayak and Chaudhuri, 2009)		$-11.64$
DF5a    (Ref. Abe et al., 2014)		$-8.80$
DF-CCSD5e    (Ref. Abe et al., 2014)		$-11.17$
NR-MRCI6    (Ref. Meyer, Bohn, and Deskevich, 2006)		$-21$
	$W_{\text{d},\text{I}}$		$W_{\text{d},\text{II}}$
GHF-ZORA (this work)	$-11.6$		$-11.4$
GKS-ZORA/B3LYP (this work)	$-10.0$		$-9.9$

Table 4. Table 4: Comparison of literature data of the 𝒫 , 𝒯 𝒫 𝒯 \mathcal{P,T} -odd eEDM interaction parameter W d subscript 𝑊 d W_{\text{d}} of the spin-rotational Hamiltonian of HgF calculated with different methods and with a quasi-relativistic GHF/GKS-ZORA approach.

Method	$\frac{W_{\text{d}}e\cdot \text{cm}}{{10}^{24}\text{Hz}\cdot h}$
Exp.+SE1   (Ref. Knight, 1981; Kozlov, 1985; Dmitriev et al., 1992)		$-47$
GRECP/RASSCF2   (Ref. Dmitriev et al., 1992)		$-48$
NR-MRCI3   (Ref. Meyer, Bohn, and Deskevich, 2006)		$-33$
DF4a   (Ref. Prasannaa et al., 2015)		$-50.42$
DF-CCSD4b   (Ref. Prasannaa et al., 2015)		$-55.81$
	$W_{\text{d},\text{I}}$		$W_{\text{d},\text{II}}$
GHF-ZORA (this work)	$-66.4$		$-65.1$
GKS-ZORA/B3LYP (this work)	$-51.1$		$-50.1$

Equations

$$\begin{gathered} \left[{\mathbf{\gamma }}^{\mu }\left(\text{i}\hslash {\partial }_{\mu }+\frac{e}{c}{A}_{\mu }\right)-m_{\text{e}}c{1}_{4\times 4}\right]\psi \\ =-\frac{d_{\text{e}}}{2}{\mathbf{\gamma }}^5{\underline{\underline{\mathbf{\sigma }}}}^{\mu \nu }{F}_{\mu \nu }\psi . \end{gathered}$$

$$\begin{array}{c}{\mathbf{\gamma }}^0=\left(\begin{array}{cc}{1}_{2\times 2}& 0_{2\times 2}\\ 0_{2\times 2}& -1_{2\times 2}\end{array}\right),\text{leavevmode}\text{leavevmode}{\mathbf{\gamma }}^k=\left(\begin{array}{cc}{0}_{2\times 2}& {\mathbf{\sigma }}^k\\ -{\mathbf{\sigma }}^k& 0_{2\times 2}\end{array}\right)\\ {\mathbf{\gamma }}^5=\text{i}{\mathbf{\gamma }}^0{\mathbf{\gamma }}^1{\mathbf{\gamma }}^2{\mathbf{\gamma }}^3\end{array}$$

$${\hat{H}}_{\text{eEDM}}=-d_{\text{e}}{\mathbf{\gamma }}^0\left[\vec{\mathbf{\Sigma }}\cdot \vec{\mathcal{E}}+\text{i}\vec{\mathbf{\alpha }}\cdot \vec{\mathcal{B}}\right].$$

$${\hat{H}}_{\text{eEDM}}\approx -d_{\text{e}}{\mathbf{\gamma }}^0\vec{\mathbf{\Sigma }}\cdot \vec{\mathcal{E}}.$$

$${\hat{H}}_{\text{eEDM}}=-d_{\text{e}}\vec{\mathbf{\Sigma }}\cdot \vec{\mathcal{E}}(\vec{r})-d_{\text{e}}\left({\mathbf{\gamma }}^0-1\right)\vec{\mathbf{\Sigma }}\cdot \vec{\mathcal{E}}(\vec{r}),$$

$$\left({\mathbf{\gamma }}^0-1\right)=-2\cdot \left(\begin{array}{cc}{0}_{2\times 2}& 0_{2\times 2}\\ 0_{2\times 2}& 1_{2\times 2}\end{array}\right).$$

$${\hat{H}}_{\text{tot}}=\hat{H}+{\hat{H}}_{\text{eEDM}}$$

$${\hat{H}}_{\text{eEDM}}=\left[\hat{P},\hat{H}\right]+{\hat{H}}_{\text{d}}.$$

$${\hat{P}}_{\text{I}}\equiv -\frac{\text{i}{d}_{\text{e}}}{\hslash e}\sum _{i=1}^{N_{\text{elec}}}\vec{\mathbf{\Sigma }}\cdot {\hat{\vec{p}}}_i$$

$$\Rightarrow \text{leavevmode}{\hat{H}}_{\text{d},\text{I}}=-d_{\text{e}}\sum _{i=1}^{N_{\text{elec}}}\left({\mathbf{\gamma }}^0-1\right)\vec{\mathbf{\Sigma }}\cdot \vec{\mathcal{E}}({\vec{r}}_i),$$

$${\hat{P}}_{\text{II}}\equiv -\frac{\text{i}{d}_{\text{e}}}{\hslash e}\sum _{i=1}^{N_{\text{elec}}}{\mathbf{\gamma }}^0\vec{\mathbf{\Sigma }}\cdot {\hat{\vec{p}}}_i$$

$$\Rightarrow \text{leavevmode}{\hat{H}}_{\text{d},\text{II}}=\frac{2\text{i}c{d}_{\text{e}}}{\hslash e}\sum _{i=1}^{N_{\text{elec}}}{\mathbf{\gamma }}^0{\mathbf{\gamma }}^5{\hat{\vec{p}}}_i^2.$$

$${\vec{\mathcal{E}}}_{\text{int}}({\vec{r}}_i)\equiv \frac{1}{e}\vec{\nabla }{V}_{\text{C}}({\vec{r}}_i),$$

$$\vec{\mathcal{E}}({\vec{r}}_i)=\sum _{A=1}^{N_{\text{nuc}}}\frac{Z_{A}e}{4\pi {ϵ}_0}\frac{{\vec{r}}_i-{\vec{r}}_A}{{\left|{\vec{r}}_i-{\vec{r}}_A\right|}^3}-\sum _{j\ne i}^{N_{\text{elec}}}\frac{e}{4\pi {ϵ}_0}\frac{{\vec{r}}_i-{\vec{r}}_j}{{\left|{\vec{r}}_i-{\vec{r}}_j\right|}^3},$$

$$\vec{\mathcal{E}}({\vec{r}}_i)\approx \sum _{A=1}^{N_{\text{nuc}}}\frac{Z_{A}e}{4\pi {ϵ}_0}\frac{{\vec{r}}_i-{\vec{r}}_A}{{\left|{\vec{r}}_i-{\vec{r}}_A\right|}^3}.$$

$$\left(\begin{array}{cc}\hat{V}({\vec{r}}_i)-ϵ& c\vec{\mathbf{\sigma }}\cdot {\hat{\vec{\pi }}}_i\\ c\vec{\mathbf{\sigma }}\cdot {\hat{\vec{\pi }}}_i& \hat{V}({\vec{r}}_i)-ϵ-2m_{\text{e}}{c}^2+2d_{\text{e}}\vec{\mathbf{\sigma }}\cdot \vec{\mathcal{E}}({\vec{r}}_i)\end{array}\right)\left(\begin{array}{c}\varphi ({\vec{r}}_i)\\ \chi ({\vec{r}}_i)\end{array}\right)=\left(\begin{array}{c}\vec{0}\\ \vec{0}\end{array}\right),$$

$$\begin{gathered} \chi (\vec{r})={\underset{\mathbf{M}}{\underbrace{\left(\underset{\mathbf{A}}{\underbrace{\left(2m_{\text{e}}{c}^2-\hat{V}(\vec{r})+ϵ\right)1_{2\times 2}}}\underset{{\text{d}}_{\text{e}}\mathbf{B}}{\underbrace{-2d_{\text{e}}\vec{\mathbf{\sigma }}\cdot \vec{\mathcal{E}}(\vec{r})}}\right)}}}^{-1} \\ \times c\left(\vec{\mathbf{\sigma }}\cdot \hat{\vec{\pi }}\right)\varphi (\vec{r}). \end{gathered}$$

$$(\mathbf{A}+d_{\text{e}}\mathbf{B}{)}^{-1}={\mathbf{A}}^{-1}(1_{2\times 2}+d_{\text{e}}\mathbf{B}{\mathbf{A}}^{-1}{)}^{-1}.$$

$$(1_{2\times 2}+d_{\text{e}}\mathbf{B}{\mathbf{A}}^{-1}{)}^{-1}=\sum _{m=0}^{\infty }{\left[-d_{\text{e}}\mathbf{B}{\mathbf{A}}^{-1}\right]}^m,$$

$$\begin{gathered} {\hat{H}}_{\text{tot,I}}^{\text{ESC}}=\hat{V}(\vec{r})1_{2\times 2}+c^2\left(\vec{\mathbf{\sigma }}\cdot \hat{\vec{\pi }}\right){\left(2m_{\text{e}}{c}^2-\hat{V}(\vec{r})+ϵ\right)}^{-1} \\ \times \sum _{m=0}^{\infty }{\left[2d_{\text{e}}\vec{\mathbf{\sigma }}\cdot \vec{\mathcal{E}}(\vec{r}){\left(2m_{\text{e}}{c}^2-\hat{V}(\vec{r})+ϵ\right)}^{-1}\right]}^m\left(\vec{\mathbf{\sigma }}\cdot \hat{\vec{\pi }}\right). \end{gathered}$$

$${\hat{H}}_{\text{tot,I}}^{\text{ESC}}={\hat{H}}_0^{\text{ESC}}+\underset{{\hat{H}}_{\text{d},\text{I}}^{\text{ESC}}}{\underbrace{c^2\left(\vec{\mathbf{\sigma }}\cdot \hat{\vec{\pi }}\right){\left(2m_{\text{e}}{c}^2-\hat{V}(\vec{r})+ϵ\right)}^{-1}2d_{\text{e}}\vec{\mathbf{\sigma }}\cdot \vec{\mathcal{E}}(\vec{r}){\left(2m_{\text{e}}{c}^2-\hat{V}(\vec{r})+ϵ\right)}^{-1}\left(\vec{\mathbf{\sigma }}\cdot \hat{\vec{\pi }}\right)}}.$$

$$\begin{gathered} {\left(2m_{\text{e}}{c}^2+ϵ-\hat{V}(\vec{r})\right)}^{-1} \\ =\frac{1}{2m_{\text{e}}{c}^2-\hat{V}(\vec{r})}\sum _{k=0}^{\infty }{\left[\frac{-ϵ}{2m_{\text{e}}{c}^2-\hat{V}(\vec{r})}\right]}^k. \end{gathered}$$

$$\begin{gathered} {\hat{H}}_{\text{d},\text{I}}^{\text{ZORA}}=c^2\left(\vec{\mathbf{\sigma }}\cdot \hat{\vec{\pi }}\right){\left(2m_{\text{e}}{c}^2-\hat{V}(\vec{r})\right)}^{-1} \\ \times 2d_{\text{e}}\vec{\mathbf{\sigma }}\cdot \vec{\mathcal{E}}(\vec{r}){\left(2m_{\text{e}}{c}^2-\hat{V}(\vec{r})\right)}^{-1}\left(\vec{\mathbf{\sigma }}\cdot \hat{\vec{\pi }}\right). \end{gathered}$$

$${\omega }_{\text{d},\text{I}}(\vec{r})=\frac{2d_{\text{e}}{c}^2}{{\left(2m_{\text{e}}{c}^2-\hat{V}(\vec{r})\right)}^2}$$

$${\hat{H}}_{\text{d},\text{I}}^{\text{ZORA}}=\left(\vec{\mathbf{\sigma }}\cdot \hat{\vec{\pi }}\right){\omega }_{\text{d},\text{I}}(\vec{r})\vec{\mathbf{\sigma }}\cdot \vec{\mathcal{E}}(\vec{r})\left(\vec{\mathbf{\sigma }}\cdot \hat{\vec{\pi }}\right).$$

$$H_{\text{d},\text{I},\lambda \rho }^{\text{ZORA}}$$

$$=⟨{\phi }_{\lambda }⟩\left(\vec{\mathbf{\sigma }}\cdot \hat{\vec{p}}\right){\omega }_{\text{d},\text{I}}(\vec{r})\vec{\mathbf{\sigma }}\cdot \vec{\mathcal{E}}(\vec{r})\left(\vec{\mathbf{\sigma }}\cdot \hat{\vec{p}}\right){\phi }_{\rho }.$$

$$H_{\text{d},\text{I},\lambda \rho }^{\text{ZORA}}=\int \text{d}\vec{r}{\left(\vec{\mathbf{\sigma }}\cdot {\vec{p}}_{\lambda }\right)}^{\ast }{\omega }_{\text{d},\text{I}}(\vec{r})\vec{\mathbf{\sigma }}\cdot \vec{\mathcal{E}}(\vec{r})\left(\vec{\mathbf{\sigma }}\cdot {\vec{p}}_{\rho }\right),$$

$$\left(\vec{\mathbf{\sigma }}\cdot \vec{v}\right)\left(\vec{\mathbf{\sigma }}\cdot \vec{u}\right)=\vec{v}\cdot \vec{u}{1}_{2\times 2}+\text{i}\vec{\mathbf{\sigma }}\cdot \left(\vec{v}\times \vec{u}\right)$$

$$\begin{gathered} H_{\text{d},\text{I},\lambda \rho }^{\text{ZORA}}=\int \text{d}\vec{r}[\text{i}{\left({\vec{p}}_{\lambda }\right)}^{\ast }{\omega }_{\text{d},\text{I}}(\vec{r})\cdot \left(\vec{\mathcal{E}}(\vec{r})\times {\vec{p}}_{\rho }\right)1_{2\times 2} \\ +{\left(\vec{\mathbf{\sigma }}\cdot {\vec{p}}_{\lambda }\right)}^{\ast }{\omega }_{\text{d},\text{I}}(\vec{r})\vec{\mathcal{E}}(\vec{r})\cdot {\vec{p}}_{\rho } \\ -\vec{\mathbf{\sigma }}\cdot \left({\left({\vec{p}}_{\lambda }\right)}^{\ast }{\omega }_{\text{d},\text{I}}(\vec{r})\times \left(\vec{\mathcal{E}}(\vec{r})\times {\vec{p}}_{\rho }\right)\right)]. \end{gathered}$$

$$\begin{gathered} H_{\text{d},\text{I},\lambda \rho }^{\text{ZORA}}={\hslash }^2\int \text{d}\vec{r}[\text{i}{\vec{\nabla }}_{\lambda }^{\ast }{\omega }_{\text{d},\text{I}}(\vec{r})\cdot \left(\vec{\mathcal{E}}(\vec{r})\times {\vec{\nabla }}_{\rho }\right) \\ +\left(\vec{\mathbf{\sigma }}\cdot {\vec{\nabla }}_{\lambda }^{\ast }\right){\omega }_{\text{d},\text{I}}(\vec{r})\left(\vec{\mathcal{E}}(\vec{r})\cdot {\vec{\nabla }}_{\rho }\right) \\ -{\vec{\nabla }}_{\lambda }^{\ast }{\omega }_{\text{d},\text{I}}(\vec{r})\left(\vec{\mathbf{\sigma }}\cdot \vec{\mathcal{E}}(\vec{r})\right){\vec{\nabla }}_{\rho } \\ +\left({\vec{\nabla }}_{\lambda }^{\ast }\cdot \vec{\mathcal{E}}(\vec{r})\right){\omega }_{\text{d},\text{I}}(\vec{r})\left(\vec{\mathbf{\sigma }}\cdot {\vec{\nabla }}_{\rho }\right)], \end{gathered}$$

$${\vec{\nabla }}_{\lambda }=\vec{\nabla }{\phi }_{\lambda };\text{leavevmode}\text{leavevmode}\text{leavevmode}{\vec{\nabla }}_{\lambda }^{\ast }=\vec{\nabla }{\phi }_{\lambda }^{\ast }$$

Full text

I Introduction

Violations of fundamental symmetries, such as those related to a combined charge conjugation ($\mathcal{C}$) and parity ($\mathcal{P}$) operation ($\mathcal{CP}$-violation), provide stringent tests for physics beyond the Standard Model of particle physics, which is often referred to as new physics.Gross (1996); Fortson, Sandars, and Barr (2003) The permanent electric dipole moment (EDM) of particles, which is the target of the present work, emerges from violation of both parity ($\mathcal{P}$) and time-reversal ($\mathcal{T}$) symmetry that is related to $\mathcal{CP}$-violation via the $\mathcal{CPT}$-theoremSchwinger (1951); Pauli, Weisskopf, and Rosenfeld (1955); Lüders (1954).Khriplovich and Lamoreaux (1997)

Even before the first experimental evidence of $\mathcal{CP}$-violation in Kaon decays,Christenson et al. (1964) Salpeter studied the effect of a permanent EDM of an electron (eEDM) and presented first calculations for hydrogen-like atoms in 1958.Salpeter (1958) Sanders found strong relativistic enhancement of effects due to an eEDM in the late 1960sSandars (1965, 1966, 1968a, 1968b) and suggested the use of polar diatomic molecules that contain a high $Z$ element (with $Z$ being the nuclear charge) for the search of the proton EDM.Sandars (1967) In the 1970s Labzowsky,Labzowsky (1978) Gorshkov et. alGorshkov, Labzowsky, and Moskalev (1982) and Sushkov and FlambaumSushkov and Flambaum (1978) showed that effects become large in polar heavy diatomic radicals. Some years later Sushkov, Flambaum and KhriplovichSushkov, Flambaum, and Khriplovich (1984); Flambaum and Khriplovich (1985) suggested to exploit small $\Omega$-doubling effects in ${}^{2}\Sigma_{1/2}$-ground states of molecules such as HgF or BaF, which were first studied theoretically by Kozlov in a semi-empirical model in 1985.Kozlov (1985) Recently, in 2014, the currently lowest upper limit on the eEDM ($d_{\text{e}}\lesssim 8.7\times 10^{-29}\leavevmode\nobreak\ e\cdot$ cm) was measured using the $H^{3}\Delta_{1}$-state of ThO.Baron et al. (2014) As $\mathcal{CP}$-violation in the Standard Model is only embedded on the level of quarks and gives rise to an EDM in the lepton sector only via higher order radiative corrections, an eEDM is a sensitive probe for new physics.Bernreuther and Suzuki (1991); Nir (2000)

The theoretical search for still more favourable candidates for eEDM experiments, employing relativistic quantum chemistry, is vivid.Kudashov et al. (2014); Skripnikov and Titov (2015); Prasannaa et al. (2015); Denis et al. (2015); Fleig, Nayak, and Kozlov (2016); Sasmal et al. (2016a, b); Sunaga et al. (2016) So far, besides the strong $Z$-dependence of $\mathcal{CP}$-violating effects, there appears to be no thorough understanding of the mechanisms that make eEDM enhancement in molecular systems large, but first attempts in this direction have been reported recentlySunaga et al. (2017). Thus, a systematic study of $\mathcal{P,T}$-odd effects in molecules is of avail. Yet relativistic many-electron calculations are computationally demanding. The main effort in relativistic calculations stems from an explicit consideration of the small component of the Dirac wave function.Dyall and Fægri, Jr. (2007) Quasi-relativistic approaches, such as the zeroth order regular approximation (ZORA), improve on computation time considerably and perform very well in most molecular calculations.van Leeuwen et al. (1994) Even in the description of relativistic molecular properties such as parity violating effects in molecules ZORA proved to be very reliable.Berger, Langermann, and van Wüllen (2005); Berger and van Wüllen (2005); Berger and Stuber (2007); Berger (2008); Nahrwold and Berger (2009); Isaev and Berger (2012, 2013, 2014)

II Methodology

II.1 Theory of eEDM Interactions in Molecules

Salpeter introduced a perturbation of the Dirac equation due to interactions with an electric dipole moment of an electron to describe permanent EDMs of atoms.Salpeter (1958) The Lorentz- and gauge-invariant formulation appears as

[TABLE]

Here $c$ is the speed of light in vacuum, $e$ is the electric constant, $m_{\text{e}}$ is the mass of the electron, $\i=\sqrt{-1}$ is the imaginary unit, $d_{\text{e}}$ is the eEDM, $\psi$ is a Dirac four-spinor, $A_{\mu}=\left(\Phi,-\vec{A}\right)$ is the four-potential with the scalar and vector potentials $\Phi$ and $\vec{A}$, $\partial_{\mu}=\frac{\partial}{\partial\mathrm{x}_{\mu}}$ is a first derivative of the four vector $\mathrm{x}_{\mu}=(ct,x,y,z)$, $\bm{\gamma}^{\mu}$, $\bm{\gamma}^{5}$ are the Dirac matrices in standard notation:

[TABLE]

and $\underline{\underline{\bm{\sigma}}}^{\mu\nu}=\left[\bm{\gamma}^{\mu},\bm{\gamma}^{\nu}\right]_{-}/2$, with $\left[a,b\right]_{-}=ab-ba$ being the commutator. $F_{\mu\nu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}$ is the field strength tensor. In the above equations index notation ($\mu,\nu=0,1,2,3;\leavevmode\nobreak\ k=1,2,3$) and Einstein’s sum convention are employed. After evaluation of the tensor product $\underline{\underline{\bm{\sigma}}}^{\mu\nu}F_{\mu\nu}$ the left-hand side of Eq. (1) reduces to a term proportional to the four-component analogue of the Pauli spin-matrix vector $\vec{\bm{\Sigma}}=\bm{1}_{2\times 2}\otimes\vec{\bm{\sigma}}$, representing the interaction with the electric field $\vec{\mathcal{E}}$ and a term proportional to $\i\vec{\bm{\alpha}}$, where $\bm{\alpha}^{k}=\bm{\gamma}^{0}\bm{\gamma}^{k}$ with $k=1,2,3$, representing the interaction with the magnetic field $\vec{{B}}$. Thus the Lorentz invariant eEDM Hamiltonian has the form

[TABLE]

It was shown that the magnetic term gives minor contributions in many-body calculations and in leading order this Hamiltonian further reduces to the interaction with the electric field:Lindroth, Lynn, and Sandars (1989)

[TABLE]

In 1963 Schiff stated, that in the non-relativistic limit the expectation value of $\hat{H}_{\text{eEDM}}$ of an atom vanishes, independent whether the elementary particles in the atom have an EDM or not.Schiff (1963) Therefore an atom in the non-relativistic limit has always a zero EDM. Addition of $0=d_{\text{e}}\left(\vec{\bm{\Sigma}}\cdot\vec{\mathcal{E}}(\vec{r})-\vec{\bm{\Sigma}}\cdot\vec{\mathcal{E}}(\vec{r})\right)$ results in the alternative formulationSandars (1965, 1966, 1968a, 1968b)

[TABLE]

with the total electrical field $\vec{\mathcal{E}}=\vec{\mathcal{E}}_{\text{int}}+\vec{\mathcal{E}}_{\text{ext}}$ being a sum of the internal and external fields $\vec{\mathcal{E}}_{\text{int}}$, $\vec{\mathcal{E}}_{\text{ext}}$. Whereas the first term on the right is zero due to Schiff’s theorem, the second appears only in the small component of the Dirac equation, because

[TABLE]

Since the momentum operator commutes with the Dirac-Hamiltonian, the effective Hamiltonian for the eEDM can be reformulated by commuting the unperturbed Hamiltonian with a modified momentum operator and defining an effective operator $\hat{H}_{\text{d}}$:

[TABLE]

In consistency with the notation of Lindroth, Lynn and SandarsLindroth, Lynn, and Sandars (1989) the modified momentum operator of the above derivation is called $\hat{P}_{\text{I}}$ and the reformulation of $\hat{H}_{\text{eEDM}}$ using $\hat{P}_{\text{I}}$ is denoted as Stratagem I, which can be summarized as

[TABLE]

Here $\hat{\vec{p}}$ is the linear momentum operator, $\hbar=h/(2\pi)$ is the reduced Planck constant and the sum runs over all $N_{\text{elec}}$ electrons of system. Additionally, introducing a factor of $\bm{\gamma}^{0}$ in the modified momentum, an alternative expression, called Stratagem II, can be derived:Mårtensson-Pendrill and Öster (1987); Commins and DeMille (2010)

[TABLE]

This operator has, within the Dirac-Coulomb picture, the advantage of being a single-particle operator, whereas $\hat{H}_{\text{d},\text{I}}$ is due to the internal electrical field, which is a function of the Coulomb potential $V_{\text{C}}$, a many-body operator (see below). Lindroth, Lynn and Sandars have pointed out, that effects of the Breit interaction and one-photon exchange, which are on the same order, namely $\mathcal{O}(\alpha^{2})$, give only minor contributions of less than one percent.Lindroth, Lynn, and Sandars (1989)

For the evaluation of $\hat{H}_{\text{d},\text{I}}$ the mathematical form of the electrical field is of interest. Starting from its definition

[TABLE]

using the molecular or atomic Coulomb potential and assuming a spherically symmetric charge distribution of the nucleus, e.g. a Gaussian charge distribution, the Gauß law is valid and the internal electrical field is

[TABLE]

with the vacuum permittivity $\epsilon_{0}$, the number of protons $Z_{A}$ of nucleus $A$, the vector in position space $\vec{r}$ and the sums running over all $N_{\text{nuc}}$ nuclei of the molecule. The first term on the right-hand side arises from the electrostatic fields of the nuclei, experienced by the electron, and the second term arises from the electrostatic fields of the other electrons. Due to the latter, $\hat{H}_{\text{d},\text{I}}$ is a many-body Hamiltonian and therefore much more difficult to treat in numerical calculation.

Fortunately, it has been shown, that the two-electron contribution is on the order of one percent only in relativistic many-body calculations and thus it is typically well justified to drop this termLindroth, Lynn, and Sandars (1989). The internal electrical field reduces then to the nuclear contribution only:

[TABLE]

In this approximation $\hat{H}_{\text{d},\text{I}}$ is a single-particle operator and for this reason there is no longer an advantage in using $\hat{H}_{\text{d},\text{II}}$ instead, even within the Dirac-Coulomb picture.

In this work both forms will be used to derive a quasi-relativistic theory in the zeroth order regular approximation (ZORA) framework. This does not only allow a comparison of the two transformations of the perturbed Dirac equation, but also provides a test of the newly developed quasi-relativistic approach, as both stratagems should yield approximately the same results.

II.2 Derivation of the eEDM ZORA equations

In the following we will derive the ZORA eEDM interaction Hamiltonian, starting from the Stratagem I Hamiltonian $\hat{H}_{\text{d},\text{I}}$ as perturbation to the molecular Dirac equation (see Equation (10)) and repeat the derivation afterwards for $\hat{H}_{\text{d},\text{II}}$, received from Stratagem II (see Equation (12)).

II.2.1 Derivation starting from Stratagem I

The $\hat{H}_{\text{d},\text{I}}$ perturbed molecular Dirac equation can be written in block matrix form as

[TABLE]

where $\phi$, $\chi$ are the large and small components of the Dirac spinor with energy $\epsilon$ shifted by $m_{\text{e}}c^{2}$, $\hat{\vec{\pi}}=\hat{\vec{p}}+e\vec{A}$ is the minimal coupling relation for the electron with $\vec{A}$ being the vector potential, $\hat{V}$ is the scalar potential energy operator and $i$ is the index of the electron, which will be dropped in the following for better readability.

Now using the elimination of small component (ESC) method,Chang, Pelissier, and Durand (1986) the small component can be expressed as

[TABLE]

This equation is valid only if the matrix $\bf M$ is invertible. This is commonly considered to be the case, because $\bf B$ has the form of the Pauli matrices (which are invertible) and $\bf A$ is a diagonal matrix and therefore does not change the invertibility of $\bf B$ by addition provided that the parameters of the fields involved are assumed to be confined appropriately.

The inversion of the matrix $\bf M$ can now be rewritten in several ways. The result of the inversion should have the following properties: (i) it should be divisible into an unperturbed and a perturbed Hamiltonian, (ii) it should be expandable in $d_{\text{e}}$ and (iii) the perturbing leading order term should be linear in $d_{\text{e}}$.

This can be achieved by first extracting ${\bf A}^{-1}=\left(2m_{\text{e}}c^{2}-\hat{V}(\vec{r})-\epsilon\right)^{-1}\bm{1}_{2\times 2}$, using the matrix relation

[TABLE]

Now the inverse expression can be expanded in $d_{\text{e}}$ when for typical field situations $\left\lvert d_{\text{e}}{\bf B}{\bf A}^{-1}\right\rvert<\bm{1}_{2\times 2}$ is assumed:

[TABLE]

and we can write down an expression for the ESC Hamiltonian of infinite order in $d_{\text{e}}$:

[TABLE]

This Hamiltonian can be separated into an unperturbed ESC-Hamiltonian $\hat{H}^{\text{ESC}}_{0}$ ($m=0$ term and potential energy term) and a perturbation due to the eEDM. The latter can be reduced to the term linear in $d_{\text{e}}$, as $\left\lvert d_{\text{e}}\right\rvert$ is very small.

[TABLE]

For regular approximationChang, Pelissier, and Durand (1986); van Lenthe, Baerends, and Snijders (1993) we extract $\left(2m_{\text{e}}c^{2}-\hat{V}(\vec{r})\right)^{-1}$ from $\left(2m_{\text{e}}c^{2}-\hat{V}(\vec{r})-\epsilon\right)^{-1}$ and expand in the orbital energy $\epsilon$:

[TABLE]

The ZORA eEDM Hamiltonian, linear in $d_{\text{e}}$, reads

[TABLE]

As $\left[\left(2m_{\text{e}}c^{2}-\hat{V}(\vec{r})\right)^{-1},\vec{\bm{\sigma}}\cdot\vec{\mathcal{E}}(\vec{r})\right]_{-}=0$, a modified ZORA factor can be defined as

[TABLE]

and the final expression for the ZORA eEDM interaction Hamiltonian reads

[TABLE]

This operator is approximately a one-electron operator (see Eq. (15)). From here on we drop now terms depending on the vector potential (thus assuming no magnetic interactions). Then the matrix elements within a one-electron basis set $\{\varphi_{\lambda}\}$ are of the form

[TABLE]

The momentum operators in position space are differential operators and therefore make simplifications complicated. These operators are hermitian and using this property, we let the left operator act on the left basis function and the right on the right basis function to receive the integral

[TABLE]

where the momentum operator acting on the basis function was simplified as $\hat{\vec{p}}\varphi_{\lambda}=\vec{p}_{\lambda}$, which is no longer an operator. Hence the commutators of all appearing elements in the above product are zero: $\left[\left(\vec{\bm{\sigma}}\cdot\vec{p}_{\rho}\right),\left(\vec{\bm{\sigma}}\cdot\vec{\mathcal{E}}\right)\right]_{-}=\left[\omega_{\text{d},\text{I}}(\vec{r}),\left(\vec{\bm{\sigma}}\cdot\vec{\mathcal{E}}\right)\right]_{-}=\left[\omega_{\text{d},\text{I}}(\vec{r}),\left(\vec{\bm{\sigma}}\cdot\vec{p}_{\rho}\right)\right]_{-}=0$. Using the Dirac relation

[TABLE]

twice (for spin-independent $\vec{u}$ and $\vec{v}$), the interaction matrix element is divided in three parts, that is

[TABLE]

This expression can be further simplified using the Graßmann identity $\vec{a}\times(\vec{b}\times\vec{c})=\vec{b}(\vec{a}\cdot\vec{c})-\vec{c}(\vec{a}\cdot\vec{b})$, receiving

[TABLE]

where the momentum operator $\hat{\vec{p}}=-\i\hbar\vec{\nabla}$ was inserted and the notation

[TABLE]

was introduced. Now we can separate terms with respect to the spatial components of the spin and write the eEDM Hamiltonian in the form

[TABLE]

This results in spin-free matrix-elements of the ZORA eEDM Hamiltonian

[TABLE]

and matrix-elements corresponding to the three spatial directions of spin $k\not=m\not=l\leavevmode\nobreak\ \wedge\leavevmode\nobreak\ k,m,l\in\left\{1,2,3\right\}$

[TABLE]

where the notation

[TABLE]

was introduced. The internal electrical field is calculated in the approximation of Eq. (15) and in the modified eEDM ZORA factor a model potential, introduced by van Wüllen,van Wüllen (1998) is used to alleviate the gauge dependence of ZORA. With that the matrix elements are implemented as

[TABLE]

[TABLE]

These integrals can not be solved analytically and therefore numerical integration on a Becke grid is used for calculations of the eEDM matrix elements. Finally the eEDM interaction can be calculated using the corresponding density and spin-density matrices ${\bf D}^{(\mu)}$, which are obtained from the molecular orbital coefficients $C_{\lambda}^{(\alpha)},C_{\lambda}^{(\beta)}$ of orbital $\lambda$ for up- ($\alpha$) and down-spin ($\beta$) by the formulas

[TABLE]

where the sum runs over all $N_{\text{occ}}$ occupied orbitals $i$. Then the molecular expectation value of the eEDM interaction Hamiltonian reads

[TABLE]

These formulas are the working equations used to calculate the expectation value of ${H}^{\text{ZORA}}_{\text{d},\text{I}}$ in first order perturbation theory. In the present calculations the (spin-)density matrices are obtained from a self-consistent field (SCF) calculation within the complex generalized Hartree-Fock (cGHF) or complex generalized Kohn-Sham (cGKS) formalism using the ZORA Hamiltonian, but could in principle be also obtained from a more sophisticated electronic structure treatment.

II.2.2 Derivation starting from Stratagem II

With the perturbation $\hat{H}_{\text{d},\text{II}}$ the block matrix form of the molecular Dirac equation looks like

[TABLE]

Again in the following the electronic index will be dropped. ESC yields for the small component

[TABLE]

Now the perturbation only appears in the numerator and no special considerations respective to the denominator (in parentheses) are needed. Hence the ESC Hamiltonian can be written down as

[TABLE]

Expanding the second term on the right one obtains an unperturbed ESC Hamiltonian, a perturbations linear in $d_{\text{e}}$ and a quadratic term in $d_{\text{e}}$:

[TABLE]

When considering the terms in first order of $d_{\text{e}}$ only and carrying out the ZORA expansion, one obtains the resulting eEDM interaction Hamiltonian as (again with the vector potential depending terms dropped)

[TABLE]

where the modified ZORA factor is defined as

[TABLE]

Its matrix elements in an arbitrary one-electron basis $\{\varphi_{\lambda}\}$ are

[TABLE]

Due to the simple form of the operator we can directly proceed in writing the matrix elements of $\hat{H}^{\text{ZORA}}_{\text{d},\text{II}}$ in terms of Eq. (32) and get in analogy to (37)

[TABLE]

and there is no spin-free part in the Hamiltonian. Again the integrals are evaluated numerically due to the appearance of the model potential $\tilde{V}(\vec{r})$ in the denominator. Finally, the total ZORA eEDM interaction energy derived from Stratagem II is evaluated from the (spin-)density matrices as presented in Eq. (39) for $\hat{H}^{\text{ZORA}}_{\text{d},\text{I}}$.

III Computational Details

For the calculation of two-component wave functions at the GHF-/GKS-ZORA level a modified versionBerger and van Wüllen (2005); Berger, Langermann, and van Wüllen (2005); Isaev and Berger (2012); Nahrwold and Berger (2009) of the quantum chemistry program package TurbomoleAhlrichs et al. (1989) was used. In order to calculate the $\mathcal{P,T}$-odd eEDM interaction, the program was extended with the ZORA eEDM Hamiltonians implemented as shown in Equations (36), (37), (39) and (47).

For density functional theory (DFT) calculations within the Kohn-Sham framework the hybrid Becke three parameter exchange functional and Lee, Yang and Parr correlation functional (B3LYP)Stephens et al. (1994); Vosko, Wilk, and Nuisar (1980); Becke (1988); Lee, Yang, and Parr (1988) was employed. For all calculations an atom centered basis set of 37 s, 34 p, 14 d and 9 f uncontracted Gaussian functions with the exponential coefficients $\alpha_{i}$ composed as an even-tempered series as $\alpha_{i}=a\cdot b^{N-i};\leavevmode\nobreak\ i=1,\dots,N$, with $b=2$ for s- and p-function and with $b=(5/2)^{1/25}\times 10^{2/5}\approx 2.6$ for d- and f-functions was used for the heavy atom (for details see Supplementary Material). This basis has been proven successful in calculations of nuclear-spin dependent $\mathcal{P}$-violating interactions in heavy polar diatomic molecules.Isaev and Berger (2012, 2013, 2014) The basis set centered at the fluorine atom was represented by an uncontracted atomic natural orbital (ANO) basis of triple-$\zeta$ qualityRoos et al. (2004). The ZORA-model potential $\tilde{V}(\vec{r})$ was employed with additional dampingLiu et al. (2002) as proposed by van Wüllen.van Wüllen (1998)

For the calculations of two-component wave functions and properties a finite nucleus was used, described by a spherical Gaussian charge distribution $\rho_{\alpha}(\vec{r})=\rho_{0}\mathrm{e}^{-\frac{3}{2\zeta_{\alpha}}\vec{r}^{2}}$, where $\rho_{0}=eZ\left(\frac{3}{2\pi\zeta_{\alpha}}\right)^{3/2}$ and the root mean square radius $\zeta_{\alpha}$ of nucleus $\alpha$ was used as suggested by Visscher and Dyall.Visscher and Dyall (1997) The mass numbers $A$ were chosen as nearest integer to the standard relative atomic mass, i.e. 19F, 137Ba, 173Yb, 201Hg, 226Ra.

The nuclear distances were optimized at the levels of GHF-ZORA and GKS-ZORA/B3LYP, respectively. For structure optimizations at the DFT level the nucleus was approximated as a point charge. The distances obtained are given in the results section.

IV Results and Discussion

In the following we focus on select diatomic molecules with a ${}^{2}\Sigma_{1/2}$-ground state, which have been studied extensively in literatureKozlov (1985); Dmitriev et al. (1992); Kozlov and Ezhov (1994); Kozlov and Labzowsky (1995); Kozlov et al. (1997); Titov, Mosyagin, and Ezhov (1996); Mosyagin, Kozlov, and Titov (1998); Quiney, Skaane, and Grant (1998); Parpia (1998); Kozlov and Derevianko (2006); Nayak and Chaudhuri (2006a, b, 2007, 2009); Kudashov et al. (2014); Abe et al. (2014); Sunaga et al. (2016); Meyer, Bohn, and Deskevich (2006); Prasannaa et al. (2015); Sasmal et al. (2016a, b), namely BaF, YbF, HgF and RaF. These systems are well suited for the validation of the here presented ZORA approach.

The eEDM contribution to the effective spin-rotational Hamiltonian for diatomic molecules in a ${}^{2}\Sigma_{1/2}$-state has the formDmitriev et al. (1992); Kozlov and Labzowsky (1995)

[TABLE]

where $\Omega=\vec{J_{\text{e}}}\cdot\vec{\lambda}$ is the projection of the total angular momentum of the electron $\vec{J}_{\text{e}}$ on the molecular axis, defined by the unit vector $\vec{\lambda}$ pointing from the heavy to the light nucleus and

[TABLE]

with the ZORA wave function $\tilde{\psi}$. In some publications, which will be referred to in the following, instead of $W_{\text{d}}$ only the effective electrical field $\mathcal{E}_{\text{eff}}=W_{\text{d}}\Omega$, was reported. In the tables below, we have converted these values for comparison then to $W_{\text{d}}$ using $\Omega=1/2$.

The GKS inter-nuclear distances are: 3.76 $a_{0}$ for YbF, 3.33 $a_{0}$ for HgH, 3.91 $a_{0}$ for HgF, 4.11 $a_{0}$ for BaF and 4.26 $a_{0}$ for RaF; and $\Omega$ was in the GKS framework: 0.473 for YbF, 0.497 for HgF and 0.500 for BaF and RaF. The GHF inter-nuclear distances are: 3.90 $a_{0}$ for YbF, 3.30 $a_{0}$ for HgH, 3.82 $a_{0}$ for HgF, 4.16 $a_{0}$ for BaF and 4.30 $a_{0}$ for RaF; and $\Omega$ was in the GHF framework: 0.498 for HgF and 0.500 for YbF,BaF and RaF.

IV.1 eEDM enhancement in the ${}^{2}\Sigma_{1/2}$-ground states of BaF and RaF

We start our discussion with BaF, which was studied already in the beginnings of the search for molecular $\mathcal{P,T}$-violation.Kozlov (1985) A number of calculations of $W_{\text{d}}$ or $W_{\text{d}}\Omega$ to compare to exists for this open-shell diatomic molecule in the literature.

From Table 1 on page 1 we see that the results calculated with the two different ZORA operators are equal within the given precision. This suggests that differences in the transformations made to obtain $H_{\text{d}}$ are of only minor importance within ZORA.

Early calculations with the generalized relativistic effective core potential (GRECP) method without effective operator (EO) based perturbative corrections were significantly lower in magnitude than Kozlov’s semi-empirical estimates. Changes upon inclusion of EO based corrections imply that the inclusion of spin-polarization is crucial for a good description of $\mathcal{P,T}$-odd properties.

The most recent four-component Dirac–(Hartree)–Fock (DF) based restricted active space (RAS) configuration interaction (DF-RASCI) calculations are in a very good agreement with GRECP calculations at the highest level of theory (RASSCF/EO, with RASSCF meaning restricted active space self-consistent field) and are also in good agreement with the semi-empirical results.Kozlov and Labzowsky (1995); Kozlov et al. (1997); Nayak and Chaudhuri (2006b)

The GHF- and GKS-ZORA calculations of this work are rather compared to electron-correlation calculations than (paired) DF, since the complex GHF/GKS approach already includes spin-polarization effects. The GHF results are in a very good agreement ($\sim 6\%$) with previous calculations, whereas the GKS-approach underestimates $W_{\text{d}}$($\sim 18\%$) in magnitude. Yet, comparison of GRECP/SCF/EO calculations and GRECP/RASSCF/EO, which both include spin-polarization, shows that spin-polarization effects are partially cancelled by other electron-correlation effects. This trend can also be seen by comparing GKS with GHF calculations. Thus ZORA performs amazingly well in the calculation of a purely relativistic property such as $W_{\text{d}}$, although not explicitly considering the small component.

Recent publications called attention to RaF as a promising candidate for the first measurement of $\mathcal{P}$- and $\mathcal{P,T}$-odd effects in molecules.Isaev and Berger (2013); Kudashov et al. (2014)

For RaF both methods, GHF and GKS, give results that are in good agreement with those of three types of relativistic electron-correlation calculations (both below 10% deviation, see Table 2 on page 2) with the DFT results agreeing slightly better with the literature values (below 5% error). Electron correlation effects as described on the DFT level display the same trends as observed for BaF.

IV.2 eEDM enhancement in the ${}^{2}\Sigma_{1/2}$-ground state of YbF

YbF is probably the best studied molecule with respect to molecular $\mathcal{CP}$-violation.

The values of $W_{\text{d}}$, calculated via Eq. (10) and via Eq. (12) are in an excellent agreement and deviations are smaller than one percent (see Table 3 on page 3). This confirms the approximate equivalence of the two used stratagems to calculate $W_{\text{d}}$ within ZORA.

As can be seen in Table 3 on page 3 there is a large discrepancy between the literature results with deviations of up to 30%. Whereas calculations with Kozlov’s semi-empirical model predict rather large values for the magnitude of $W_{\text{d}}$ (about $12.6\times 10^{24}\leavevmode\nobreak\ \frac{\text{Hz}\cdot h}{e\cdot\text{cm}}$),Kozlov and Ezhov (1994); Kozlov et al. (1997) early GRECP/RASSCF calculations and Dirac–Fock calculations without consideration of electron-correlation yield much lower absolute values.Titov, Mosyagin, and Ezhov (1996); Nayak and Chaudhuri (2006a); Abe et al. (2014) This may result from the neglect of spin-polarization effects, which play a major role, as has been discussed in the previous section.

Yet, more recent four-component electron-correlation calculations and GRECP calculations with perturbative effective operator corrections, which include spin-polarization, show a better agreement and can be taken as the most reliable of the shown literature values. Mosyagin, Kozlov, and Titov (1998); Nayak and Chaudhuri (2007); Nayak and Das (2009); Abe et al. (2014); Sunaga et al. (2016)

In comparison to DF, concerning the parameter $W_{\text{d}}$ the deviation of the here presented ZORA approach is below 5% for GHF and of the order of 10% for GKS.

Restricted DF calculationsQuiney, Skaane, and Grant (1998) deviate almost by a factor of two from these results and therefore it can be assumed that a different definition of the effective Hamiltonian was used. Calculations within a non-relativistic framework, reported in Ref. Meyer, Bohn, and Deskevich, 2006 overestimate the magnitude of $W_{\text{d}}$ by about a factor of two and are not reliable (see also Table 4 on page 4).

Again, electron correlation corrections as estimated on the DFT level lower the absolute value of $W_{\text{d}}$.

IV.3 eEDM enhancement in the ${}^{2}\Sigma_{1/2}$-ground state of HgF

Although not studied in such detail as YbF, there is a good amount of literature on HgF as well, partially dating back to the 1980s.Kozlov (1985) Furthermore the calculations, which will be discussed in the following, show that mercury compounds can provide large enhancements of $\mathcal{P}$- and $\mathcal{CP}$-violation and therefore are very interesting to study.

As for the other compounds discussed the agreement between $W_{\text{d},\text{I}}$ and $W_{\text{d},\text{II}}$ is excellent, confirming the validity of the two transformations of the eEDM Hamiltonian in the ZORA picture.

Whereas GRECP/RASSCF and DF results are in line with the semi-empirical estimates by Kozlov, although spin-polarization effects are not accounted for Dmitriev et al. (1992); Prasannaa et al. (2015), more recent four-component relativistic coupled cluster calculations predict about 10% larger absolute values for $W_{\text{d}}$.

As opposed to the trends observed for the molecules discussed before, for HgF the GHF-ZORA approach overestimates the magnitude of the results from four-component relativistic electron-correlation calculations by more than 15%. This appears to be caused by a very pronounced energetic splitting in the Kramers pair structure of the valence orbitals below the singly occupied orbital, which have $\sigma$-symmetry. This splitting is much smaller in YbF, RaF and BaF.

The GKS-ZORA results instead appear to be much closer to the literature data (about 8% deviation from DF-CCSD and even less from DF or GRECP/RASSCF). Here the additional electron-correlation effects as estimated on the DFT level, which lead to a reduction of the absolute value of $W_{\text{d}}$ in comparison to GHF, play a much more important role than for YbF within the GHF/GKS-ZORA approach and the difference between the GHF and GKS results is much larger for HgF.

In four-component calculations the electron correlation effects seem to be less pronounced, although accounting for spin-polarization leads to very different results when one compares results of GRECP/RASSCF with those of DF-CCSD. This may be caused by partial cancellation of spin-polarization with other correlation effects.

V Conclusion

In this paper we derived a ZORA-based perturbation Hamiltonian for the description of $\mathcal{P,T}$-odd interactions due to an electron electric dipole moment in molecules. With calculations of promising candidates for a search of an eEDM, we could show that a quasi-relativistic ZORA approach is well suited for the calculation of a purely relativistic effect, although the small component of the wave function is not considered explicitly. The accuracy for prediction of $W_{\mathrm{d}}$ is estimated to be on the order of about 20 % for the ${}^{2}\Sigma_{1/2}$-ground state molecules studied herein, if one considers recent results obtained with electron correlation approaches as a benchmark. This level of accuracy is presently fully sufficient for the identification of molecular candidate systems for an eEDM search.

With the quasi-relativistic approach presented in this work an efficient calculation of the eEDM enhancement in molecules is possible. In future work we will study a larger number of molecules with this approach in order to achieve a deeper understanding of the mechanisms that lead to sizeable $\mathcal{P,T}$-odd properties in molecules.

Acknowledgements.

We thank Timur Isaev for discussions. Financial support by the State Initiative for the Development of Scientific and Economic Excellence (LOEWE) in the LOEWE-Focus ELCH and computer time by the center for scientific computing (CSC) Frankfurt are gratefully acknowledged.