Toolbox approach for quasi-relativistic calculation of molecular properties for precision tests of fundamental physics

TL;DR

This paper introduces a flexible method for calculating relativistic molecular properties using quasi-relativistic wavefunctions, enabling efficient and accurate predictions crucial for fundamental physics tests.

Contribution

It presents a general scheme for relativistic property calculations with two-component wavefunctions, demonstrated with ZORA, applicable to symmetry-violating properties in molecules.

Findings

ZORA-based results agree well with four-component calculations.

The method effectively predicts properties for molecules like RaF and TlF.

RaF's potential in experiments searching for fundamental symmetry violations is confirmed.

Abstract

A generally applicable approach for the calculation of relativistic one-electron properties with two-component wave functions is presented. The formalism is explicitly evaluated for the example of quasi-relativistic wavefunctions obtained within the zeroth order regular approximation (ZORA). The wide applicability of the scheme is demonstrated for the calculation of parity ($\mathcal{P}$) and time-reversal ($\mathcal{T}$) symmetry violating properties, which are important for searches of physics beyond the standard model of particle physics. The quality of the ZORA results is shown exemplarily for the molecules RaF and TlF by comparison to data from four-component calculations as far as available. Finally, the applicability of RaF in experiments that search for $\mathcal{P,T}$-violation not only in the electronic but also in quark sector is demonstrated.

Tables

Table 1. Table 1: Representation of various 𝒫 , 𝒯 𝒫 𝒯 \mathcal{P,T} -odd properties in the formalism of general tensor functions.

Property	Prefactor	$\hat{s}$	$\widehat{𝘁}\circ$	$\hat{\partial }\circ$	${𝚪}^{i,j,k}$
$W_d$	$\frac{2c\hslash }{e\Omega }$	-	-	$\overrightarrow{\nabla }\cdot \overrightarrow{\nabla }\otimes$	${𝚪}^{2,0,1}$
$W_{ℳ}$	$\frac{3}{2\Omega }ce{k}_{\mathrm{em}}$	$\frac{z-z_A}{{\left|\overrightarrow{r}-{\overrightarrow{r}}_A\right|}^5}$	$-(\overrightarrow{r}-{\overrightarrow{r}}_A)\times$	-	${\overrightarrow{𝚪}}^{1,(1,2,3),0}$
$W_{𝒮}$	$\frac{2\pi }{3}$	$\delta \left({\overrightarrow{r}}_A-\overrightarrow{r}\right)$	-	$\overrightarrow{\nabla }\otimes$	${𝚪}^{0,0,0}$
$W_m$	$4\left(\frac{{\mu }_N}{A_A}+\frac{{\mu }_A}{Z_A}\right)\frac{c}{\hslash }\frac{k_{\mathrm{em}}}{k}$	${\left|\overrightarrow{r}-{\overrightarrow{r}}_A\right|}^{-3}$	$(\overrightarrow{r}-{\overrightarrow{r}}_A)\times$	$\overrightarrow{\nabla }\times$	${\overrightarrow{𝚪}}^{1,(1,2,3),1}$
$W_s$	$\frac{G_F}{\sqrt{2}\Omega }$	${\rho }_{\mathrm{nuc}}\left(\overrightarrow{r}\right)$	-	-	${𝚪}^{2,0,1}$
$W_T$	$\sqrt{2}{G}_F$	${\rho }_{\mathrm{nuc}}\left(\overrightarrow{r}\right)$	-	-	${𝚪}^{2,3,1}$

Table 2. Table 2: Nuclear spin-dependent, electron spin-independent 𝒫 , 𝒯 𝒫 𝒯 \mathcal{P,T} -odd electronic structure parameters of the diamagnetic molecule 205 TlF evaluated at the level of cGHF- and cGKS-ZORA with a large even tempered basis set with (wr) and without renormalisation (wor) of the density according to eq. ( 34 ). Comparison to literature values determined with different computational method. The value μ ( 205 Tl ) = 1.6382135 μ N \mu(^{205}\mathrm{Tl})=1.6382135\leavevmode\nobreak\ \mu_{\mathrm{N}} was used for the nuclear magnetic moment of 205 Tl. Mills  et al. ( 1993 )

Method	$W_T/h\mathrm{Hz}$	$W_m/\frac{{10}^{18}h\mathrm{Hz}}{e\mathrm{cm}}$	$W_{𝒮}/a_0^4$
cGHF-ZORA-wor	$4697$	$-4.74$	$8443$
cGHF-ZORA-wr	$4690$	$-4.72$	$8428$
cGKS-ZORA-B3LYP-wr	$3375$	$-3.10$	$5720$
DHFa (Ref. Quiney et al., 1998)	$4632$	$-4.78$	$8747$
GRECP-RCC-SDb (Ref. Petrov et al., 2002)	-	$-4.04$	$7635$
DFc (Ref. Parpia, 1997)	-	$-5.46$	$7738$

Table 3. Table 3: Electron and nuclear spin-dependent 𝒫 , 𝒯 𝒫 𝒯 \mathcal{P,T} -odd electronic structure parameters of the paramagnetic molecule 223 RaF evaluated at the level of cGHF- and cGKS-ZORA with a large even tempered basis set with (wr) and without renormalisation (wor) of the density according to eq. ( 34 ). Comparison to available literature values determined with different computational methods. The value μ ( 223 Ra ) = 0.271 μ N \mu(^{223}\mathrm{Ra})=0.271\leavevmode\nobreak\ \mu_{\mathrm{N}} was used for the nuclear magnetic moment of 223 Ra. Mills  et al. ( 1993 ) In all our calculations Ω = 0.500 Ω 0.500 \Omega=0.500 .

Method	$W_T/h\mathrm{Hz}$	$W_m/\frac{{10}^{18}h\mathrm{Hz}}{e\mathrm{cm}}$	$W_{𝒮}/a_0^4$	$W_s/h\mathrm{kHz}$	$W_d/\frac{{10}^{24}h\mathrm{Hz}}{e\mathrm{cm}}$	$W_{ℳ}/\frac{{10}^{33}h\mathrm{Hz}}{ce{\mathrm{cm}}^2}$
cGHF-ZORA-wor	$-1810$	$0.66$	$-4235$	$-152$	$-27.2$	$-1.17$
cGHF-ZORA-wr	$-1809$	$0.66$	$-4229$	$-152$	$-27.2$	$-1.17$
cGKS-ZORA-B3LYP-wr	$-1617$	$0.58$	$-3686$	$-138$	$-24.6$	$-1.03$
FS-RCCSD+${\Delta }_{\text{basis}}$+${\Delta }_{\text{triples}}$a (Ref.Kudashov et al., 2014)	-	-	$-4260$	$-139$	$-25.6$	-
DF-CCSDb (Ref. Sasmal et al., 2016)	-	-	-	$-141$	$-25.4$	-

Equations

$$\psi =\left(\begin{array}{c}{\psi }^{\text{L}}\\ {\psi }^{\text{S}}\end{array}\right)=\left(\begin{array}{c}{\psi }_{\alpha }^{\text{L}}\\ {\psi }_{\beta }^{\text{L}}\\ {\psi }_{\alpha }^{\text{S}}\\ {\psi }_{\beta }^{\text{S}}\end{array}\right).$$

$$\hat{\mathbf{O}}=\left(\begin{array}{ccc}{\hat{\mathbf{O}}}_{2\times 2}^{\text{LL}}& & {\hat{\mathbf{O}}}_{2\times 2}^{\text{LS}}\\ {\hat{\mathbf{O}}}_{2\times 2}^{\text{SL}}& & {\hat{\mathbf{O}}}_{2\times 2}^{\text{SS}}\end{array}\right)$$

$$⟨\hat{\mathbf{O}}⟩=\int {\mathrm{d}}^4\vec{x}\text{leavevmode}\left[{\Omega }^{\mathrm{LL}}(\vec{x})+{\Omega }^{\mathrm{SL}}(\vec{x})+{\Omega }^{\mathrm{LS}}(\vec{x})+{\Omega }^{\mathrm{SS}}(\vec{x})\right].$$

$${\mathbf{\Gamma }}^{i,j,k}={ı}^k{\mathbf{\sigma }}^i\otimes {\mathbf{\sigma }}^j;\text{leavevmode}i,j=0,1,2,3\text{leavevmode}\wedge \text{leavevmode}k\in {\mathbb{N}}_0,$$

$${\mathbf{\Gamma }}^{0,0,k}$$

$${\mathbf{\Gamma }}^{1,0,k}$$

$${\mathbf{\Gamma }}^{2,0,k}$$

$${\mathbf{\Gamma }}^{3,0,k}$$

$${\vec{\mathbf{\Gamma }}}^{0,(1,2,3),k}$$

$${\vec{\mathbf{\Gamma }}}^{1,(1,2,3),k}$$

$${\vec{\mathbf{\Gamma }}}^{2,(1,2,3),k}$$

$${\vec{\mathbf{\Gamma }}}^{3,(1,2,3),k}$$

$$\begin{gathered} {\Gamma }^{i,j,k}({\vec{r}}_1,{\vec{r}}_1^{\prime }) \\ =\hat{\mathcal{D}}\left[{\left(\Psi (\vec{X})\right)}^{†}{\vec{\mathbf{\Gamma }}}_1^{i,j,k}\otimes {\mathbf{\Gamma }}_2^{0,0,0}\otimes \cdots \otimes {\mathbf{\Gamma }}_{N_{\mathrm{elec}}}^{0,0,0}\Psi ({\vec{X}}^{\prime })\right]. \end{gathered}$$

$$\hat{\mathcal{D}}=N_{\mathrm{elec}}\mathfrak{Re}\int\limits_{\sigma_{1}=\sigma_{1}^{\prime}}\mathrm{d}\sigma_{1}\idotsint\limits_{\vec{x}_{i}=\vec{x}^{\prime}_{i};\leavevmode\nobreak\ i=1,\dots,N}\mathrm{d}^{4}\vec{x}_{2}\cdots\mathrm{d}^{4}\vec{x}_{N},$$

$${\Gamma }^{0,j,k}(\vec{r},{\vec{r}}^{\prime })$$

$${\Gamma }^{1,j,k}(\vec{r},{\vec{r}}^{\prime })$$

$${\Gamma }^{2,j,k}(\vec{r},{\vec{r}}^{\prime })$$

$${\Gamma }^{3,j,k}(\vec{r},{\vec{r}}^{\prime })$$

$$\mathbf{\varsigma }(j,k)={ı}^k{\mathbf{\sigma }}^j\text{leavevmode}.$$

$$\begin{gathered} {\rho }^{IJ,j,k}({\vec{r}}_1,{\vec{r}}_1^{\prime }) \\ =\hat{\mathcal{D}}\left[{\left({\Psi }^I(\vec{X})\right)}^{†}{\mathbf{\varsigma }}_1(j,k)\otimes {\mathbf{\sigma }}_2^0\otimes \cdots \otimes {\mathbf{\sigma }}_{N_{\mathrm{elec}}}^0{\Psi }^J({\vec{X}}^{\prime })\right]. \end{gathered}$$

$$\begin{gathered} {\psi }^{\mathrm{S}}(\vec{x})=\underset{=c\mathbf{\omega }(\vec{x})}{\underbrace{c{\left(2m_{\mathrm{e}}{c}^2{1}_{2\times 2}-{\hat{\mathbf{V}}}_{\mathrm{diag}}^{\mathrm{SS}}(\vec{x})+ϵ1_{2\times 2}\right)}^{-1}}} \\ \times \underset{\mathbf{\varpi }(\vec{x})}{\underbrace{\left(\vec{\mathbf{\sigma }}\cdot \hat{\vec{p}}+{\hat{\mathbf{V}}}_{\mathrm{off}}^{\mathrm{SL}}(\vec{x})\right)}}{\psi }^{\mathrm{L}}(\vec{x}). \end{gathered}$$

$$\mathbf{\varpi }(\vec{x})=\sum _{j=0}^3\sum _{k=0}^1\mathbf{\varsigma }(j,k){\hat{\pi }}^{j,k}(\vec{x})$$

$$\mathbf{\varpi }\approx \vec{\mathbf{\sigma }}\cdot \hat{\vec{p}}=-\hslash ı\vec{\mathbf{\sigma }}\cdot \vec{\nabla },$$

$${\psi }^{\text{L}}(\vec{x})$$

$${\psi }^{\text{S}}(\vec{x})$$

$${\Psi }^{\text{ZORA}}=\left|{\psi }_1^{\text{ZORA}}\cdots {\psi }_N^{\text{ZORA}}\right|$$

$${\psi }_i^{\text{ZORA}}=\sum _{\mu =1}^{N_{\text{basis}}}{\vec{C}}_{i\mu }{\chi }_{\mu },$$

$${\vec{C}}_{i\mu }=\left(\begin{array}{c}{C}_{i\mu }^{\alpha }\\ C_{i\mu }^{\beta }\end{array}\right).$$

$${\rho }^{IJ,j,k}(\vec{r},{\vec{r}}^{\prime })=\sum _{i=1}^{N_{\text{orb}}}{n}_i{\rho }_i^{IJ,j,k}(\vec{r},{\vec{r}}^{\prime }).$$

$${\rho }_i^{\mathrm{LL},j,k}(\vec{r},{\vec{r}}^{\prime })=$$

$${\rho }_i^{\mathrm{LS},j,k}(\vec{r},{\vec{r}}^{\prime })=$$

$${\rho }_i^{\mathrm{SL},j,k}(\vec{r},{\vec{r}}^{\prime })=$$

$${\rho }_i^{\mathrm{SS},j,k}(\vec{r},{\vec{r}}^{\prime })=$$

$${\rho }_i^{IJ,j,k}(\vec{r},{\vec{r}}^{\prime })=\sum _{\mu \nu }^{N_{\mathrm{basis}}}\sum _{l\in M(I)}\sum _{m\in M(J)}{\mathcal{D}}_{i\mu \nu }^{j,k,(l,m)}{\mathcal{I}}_{\mu \nu }^{(l,m)}(\vec{r},{\vec{r}}^{\prime }),$$

$${\mathcal{D}}_{i\mu \nu }^{j,k,(l,m)}$$

$$=\mathfrak{Re}\left\{{\left({\vec{C}}_{i\mu }\right)}^{†}\mathbf{\varsigma }(l,1-{\delta }_{l0})\mathbf{\varsigma }(j,k)\mathbf{\varsigma }(m,3-3{\delta }_{m0}){\vec{C}}_{i\nu }\right\}$$

$${\mathcal{I}}_{\mu \nu }^{(l,m)}(\vec{r},{\vec{r}}^{\prime })=(\hslash c{)}^{2-{\delta }_{l0}-{\delta }_{m0}}{\partial }_{\mu }^{(l)}{\omega }^{1-{\delta }_{l0}}{\omega }^{\prime 1-{\delta }_{m0}}{\partial }_{\nu }^{\prime (m)},$$

$$D_{\mu \nu }^{(k)}=\sum _{i=1}^{N_{\text{occ}}}\underset{D_{i\mu \nu }^{(k)}}{\underbrace{{\vec{C}}_{i\mu }^{†}{\mathbf{\sigma }}^k{\vec{C}}_{i\nu }}}.$$

Full text

I Introduction

In order to explore the boundaries of the standard model of particle physics and of general relativity, a wealth of experiments was proposed in the last decades (see e.g. Ref. DeMille, Doyle, and Sushkov, 2017) aiming for the detection of new physics. Molecular systems have gained meanwhile increasing importance for this research direction, which has recently been reviewed for instance in Ref. Safronova et al., 2018. Modern experiments take advantage of the rich, but also highly complex vibronic internal structure of molecules, which is why theory plays a crucial role for design and interpretation of established and future experiments.

From this search for new physics emerges also a large number of new properties of interest, of which many violate symmetries that are either conserved or only very weakly broken in established physics (see e.g. Ref. Khriplovich and Lamoreaux, 1997). Several of the proposed phenomena are considered to be favourably enhanced in heavy-elemental molecules due to relativistic effects and therefore demand on the theory side a corresponding description of the electronic wave function, including spin-orbit coupling.

An accurate treatment of spin-orbit coupling in heavy elements requires at least a two-component wave function, which is nowadays available in many quantum chemistry software packages. However, naturally relativistic properties are typically described in a four-component framework. Therefore, all operators have to be transformed from four-component to two-component pictures in order to be consistent with quasi-relativistic wave function formulations. As most programs do not include direct support for four-component treatments, this requires derivation of analytical expressions of two-component operators and commonly also a new, often tedious and error-prone implementation for every property.

In this paper we develop a general formalism for flexible calculation of arbitrary relativistic one-electron properties in a two-component framework. We introduce our formalism in terms of two-component one-electron density functions and we provide a general formulation of relativistic one-electron properties for the case of zeroth order regular approximation (ZORA) wave functions. Our generally applicable scheme is demonstrated on the example of different sources of simultaneous parity ($\mathcal{P}$) and time-reversal ($\mathcal{T}$) symmetry violation in the diatomic molecules RaF and TlF, which are promising candidates for a first measurement of $\mathcal{P,T}$-violation beyond the Standard Model.Hinds and Sandars (1980); Hunter et al. (2012); Isaev and Berger (2013); Norrgard et al. (2017) In future applications the concepts and implementation derived in the following can also be used to study conventional nuclear magnetic resonance (NMR) shielding constants and parity violating shifts in chiral moleculesNahrwold and Berger (2009); Nahrwold, Berger, and Schwerdtfeger (2014) and a wealth of further properties.

II Theory

II.1 Relativistic one-particle wave functions and

operators

Relativistic electronic one-particle wave functions have a bi-spinor structure, which can be decomposed as

[TABLE]

Here $\psi^{\text{L}}$ and $\psi^{\text{S}}$ are the spinors of the large or upper component and the small or lower component of the Dirac one-particle wave function $\psi$, respectively. $\psi_{\alpha}$ and $\psi_{\beta}$ are spin orbitals with electron spin up and down, respectively. A similar structure is also found for one-particle operators in four component theory. Operators can be decomposed in a $2\times 2$-block structure:

[TABLE]

Thus, we can deduce four kinds of one-particle densities, which appear as contribution to all expectation values of custom one-electron operators in one-particle four-component theory. These are contributions from large component-large component (LL), large component-small component (LS), small component-large component (SL) and small component-small component (SS) matrix elements. This formulation allows a decomposition of expectation values of four-component one-electron operators in a sum of modified relativistic one-electron density functions $\Omega^{IJ}(\vec{x})=\left(\psi^{I}(\vec{x})\right)^{\dagger}\hat{\mathbf{O}}^{IJ}_{2\times 2}\psi^{J}(\vec{x})$, where $IJ$ can be LL, LS, SL or SS:

[TABLE]

In the following all $\Omega^{IJ}$ are reformulated in terms of two-component density functions exploiting the Lorentz symmetry of relativistic operators. Finally, this reformulation will allow a generally applicable procedure for computation of $\Omega^{IJ}$ within an approximate two-component theory.

II.2 Lorentz group and relativistic density functions

The explicit representation of four-component one-electron operators can be reduced to a linear-combination of a special unitary basis of matrices in four-dimensional space that corresponds to the Lorentz group $SU(2)\times SU(2)$, which represents the $2\times 2$-block structure of the Dirac equation. This basis is formed by the 16 (or rather 32, allowing an imaginary phase of $\imath^{k}$ with the imaginary unit $\imath=\sqrt{-1}$ and $k\in\mathbb{N}_{0}$) Dirac matrices:

[TABLE]

where $\bm{\sigma}^{i}$ for $i=1,2,3$ are the Pauli spin matrices and $\bm{\sigma}^{0}=\bm{1}_{2\times 2}$ is the two dimensional identity matrix, which build the basis of $SU(2)$, and $\otimes$ is the Kronecker product. We include an imaginary phase $\imath^{k}$ in order to enable time-reversal symmetry violation in our formalism, as well.

The usual notation for theses matrices is

[TABLE]

Employing this structure we can define 32 Dirac one-electron density functions of electron 1 which are build as

[TABLE]

We have introduced in this equation the many-particle wave function $\Psi$, which depends on the set $\vec{X}=\left\{\vec{x}_{1},\dots,\vec{x}_{N}\right\}$ of $N_{\mathrm{elec}}$ combined spin and spatial coordinates, with $\vec{x}_{i}=\left\{\sigma_{i},\vec{r}_{i}\right\}$ for each electron $i$, and the one-electron density projection operator

[TABLE]

which integrates out all coordinates except one spatial coordinate $\vec{r}_{1}$. Coordinates with a prime are introduced to formally allow a discrimination between the left and right function that contributes to the density function, in order to enable the definition of operators acting on the left or right function, respectively.

Employing the $2\times 2$-block structure of the one-particle Dirac equation we can decompose the Dirac one-electron density functions further by introducing four kinds of two-component densities for each of the four four-component densities $\rho^{IJ,j,k}$. The two-component density functions can be classified as i) number density functions $\rho^{IJ,0,2k}$, ii) number current density functions $\rho^{IJ,0,2k+1}$, iii) spin density functions $\vec{\rho}^{IJ,(1,2,3),2k}$ and iv) spin current density functions $\vec{\rho}^{IJ,(1,2,3),2k+1}$. Thus there are in total $4\times 4=16$ two-component one-electron density functions, which can be combined linearly to give the 32 Dirac one-electron density functions:

[TABLE]

We introduce the spin tensor function in correspondence to the four Pauli matrices with a complex phase:

[TABLE]

Using this spin tensor function the two-component one-electron density functions of electron 1 are build with the spinors of the Dirac one-particle wave function containing either only the large component for the first electron $\Psi^{\mathrm{L}}$ or only small component for the first electron $\Psi^{\mathrm{S}}$ as

[TABLE]

II.3 Quasi-relativistic approximation and ZORA

In order to be able to calculate these relativistic density functions in two-component theories we need to approximate the small component wave function. The exact transformation of the large-component one-particle wave function into the small-component one-particle wave function is given by:

[TABLE]

Here $\epsilon$ is the electronic energy of the single particle state, $\hat{\vec{p}}=-\imath\hbar\vec{\nabla}$ is the one-electron linear momentum operator in position space with $\hbar=\frac{h}{2\pi}$ being the reduced Planck’s constant, $\hat{\mathbf{V}}^{\mathrm{SS}}_{\mathrm{diag}}$ is the SS-block of the potential effective one-electron operator appearing on the diagonal of the Hamiltonian, $\hat{\mathbf{V}}^{\mathrm{SL}}_{\mathrm{off}}$ is the SL-block of a potential effective one-electron operator appearing on the off-diagonal of the Hamiltonian, as e.g. a vector potential or the Breit operator, both as effective one-electron operators. $m_{\mathrm{e}}$ is the electron mass and $c$ is the speed of light in vacuum.

For any two-component method we can replace $\psi^{\mathrm{L}}$ approximately by the wave-function optimized within this specific method and use the corresponding approximate transformation matrix $\bm{\omega}$ to construct an approximate $\psi^{\mathrm{S}}$.

As $\hat{\mathbf{V}}^{\mathrm{SL}}_{\mathrm{off}}(\vec{x})$ is an effective one-electron operator, the modified momentum matrix $\bm{\varpi}$ can be decomposed into a linear combination of spin tensor functions:

[TABLE]

with the modified one-electron momentum operators $\hat{\vec{\pi}}^{j,k}$. Neglecting all off-diagonal potentials, $\bm{\varpi}$ simplifies to

[TABLE]

so that $\hat{\vec{\pi}}^{(1,2,3),1}=-\hbar\vec{\nabla}$. For all $j=0$ or $k=0$, we have $\hat{\pi}^{j,k}_{i}=0$.

In the present paper we implemented the scheme within zeroth order regular approximation (ZORA) without potentials appearing on the off-diagonal:

[TABLE]

where we assume an unperturbed scalar transformation factor $\omega$, which depends only on a model-potential $\tilde{V}$ as introduced by van Wüllen.van Wüllen (1998) The model-potential is an effective one-electron operator that does not depend on the molecular wave function and is designed to alleviate the gauge variance of ZORA. All perturbations to the diagonal potential $\hat{\mathbf{V}}_{\mathrm{diag}}$ are neglected in $\tilde{V}$.

Corrections stemming from the presence of vector potentials appear as additional terms and can be included afterwards by calculation of matrix elements of the corresponding operators, which is straightforward in the present approach.

In the mean-field approximation the ZORA $N_{\mathrm{elec}}$-electron wave function is represented as $N_{\mathrm{elec}}$ Slater determinant of one-electron wave functions $\psi_{i}^{\text{ZORA}}$, which are composed of $N_{\mathrm{basis}}$ real spatial basis functions $\chi_{\mu}$ and complex two component coefficients $\vec{C}_{i\mu}$:

[TABLE]

with

[TABLE]

where

[TABLE]

Only the coefficients depend on the spin component, so that spin can in this form easily be integrated out. Furthermore, in this approach we can define for each molecular orbital (MO) $i$ a MO-density function $\rho^{IJ,j,k}_{i}$. The total ZORA-density function is a sum of all $N_{\text{orb}}$ MO-density functions weighted by their occupation number $n_{i}$, which is 0 for unoccupied or 1 for occupied orbitals:

[TABLE]

Within ZORA the one-electron MO-density functions $\rho^{IJ,j,k}_{i}$ are evaluated explicitly as

[TABLE]

Here we used the short-hand notations: $\chi_{\mu}(\vec{r})=\chi_{\mu}$, $\chi_{\mu}(\vec{r}^{\prime})=\chi^{\prime}_{\mu}$, $\vec{\nabla}\chi_{\mu}=\vec{\partial}_{\mu}$, $\vec{\nabla}^{\prime}\chi^{\prime}_{\mu}=\vec{\partial}^{\prime}_{\mu}$, $\omega^{\text{ZORA}}(\vec{r})=\omega$ and $\omega^{\text{ZORA}}(\vec{r}^{\prime})=\omega^{\prime}$. Generalization of eqs. (20a) to (20d) in terms of $\mathcal{D}$, $\mathcal{I}$ gives

[TABLE]

where $M(I)=\begin{cases}\left\{1,2,3\right\}&\text{if}\leavevmode\nobreak\ I=\mathrm{S}\\ \left\{0\right\}&\text{if}\leavevmode\nobreak\ I=\mathrm{L}\end{cases}$ is a mapping between spin indices $l,m$ of $\bm{\mathcal{D}}^{j,k,(l,m)}_{i}$ and $\bm{\mathcal{I}}^{(l,m)}$ and the relativistic type of the two-component density function $I,J$ and

[TABLE]

with $\partial^{(i)}_{\mu}=\begin{cases}\frac{\partial\chi_{\mu}}{\partial(\vec{r})_{i}}&\text{if}\leavevmode\nobreak\ i=1,2,3\\ \chi_{\mu}&\text{if}\leavevmode\nobreak\ i=0\end{cases}$ and analogously with prime.

Here the density matrices $\bm{\mathcal{D}}^{j,k,(l,m)}=\sum\limits_{i=1}^{N_{\mathrm{occ}}}\bm{\mathcal{D}}^{j,k,(l,m)}_{i}$ can be constructed from the four (eight when separating real and imaginary parts) two-component density matrices $\mathbf{D}^{(k)}$, which are defined in terms of Pauli matrices as

[TABLE]

The explicit expressions for $\bm{\mathcal{D}}^{j,k,(k,l)}$ in terms of $\mathbf{D}^{(k)}$ can be evaluated following eq. (22) and using commutator relations of Pauli matrices.

II.4 Generic tensor function formulation of molecular

properties

With the two-component formulation of relativistic density functions introduced in the previous section we can formulate a generic one-electron tensor function which can directly be connected to an arbitrary one-electron molecular property. For this purpose we introduce, besides the relativistic $\Gamma^{i,j,k}$ or quasi-relativistic one-electron density function $\rho^{IJ,j,k}$, a general differential tensor operator $\hat{\bm{\partial}}(\vec{r}^{\prime}\vee\vec{r})$, which acts on $\vec{r}^{\prime}$ or $\vec{r}$. Within this flexible formulation we do not use turn-over rule, which can complicate calculations as pointed out in Ref. Kutzelnigg, 1988. Furthermore a general tensor operator $\hat{\bm{\mathsf{t}}}(\vec{r})$ and a general scalar operator $\hat{s}$ are defined. All these operators may be defined arbitrarily. In the following section, some explicit realisations for these operators will be discussed. A generic one-electron tensor function can now be written as

[TABLE]

or for two-component density functions as

[TABLE]

where $\circ=\otimes\leavevmode\nobreak\ ,\cdot,\leavevmode\nobreak\ \times$, with $\otimes$ being the outer product, $\cdot$ being the inner product and $\times$ being the cross product. The latter is defined only within $\mathbb{R}^{3}$. Here $j$ can be [math] or $(1,2,3)$. In the latter case $\Gamma^{i,j,k}$ is a three dimensional vector. This formulation allows the flexible construction of one-electron operators that can correspond to molecular properties. The expectation value is received by integration over space:

[TABLE]

With explicit expressions for quasi-relativistic ZORA density functions given in eqs. (21), (22) and (23). and the definitions of the relativistic density functions in eq. (8) the working equation for the computation of the expectation value of an arbitrary relativistic one-electron properties within ZORA is

[TABLE]

Here $\mathcal{M}(i)$ maps the index $i$ of the Dirac density function to a sum of two-component density functions with appropriate sign following eqs. (8).

The above expressions are valid for properties within first order perturbation theory only. In case of second order properties the corresponding density functions have to be constructed with perturbed density matrices $\tilde{\mathbf{D}}^{(k)}$ or rather transition density matrices. These can be received via common algorithms by solving the coupled perturbed Hartree-Fock (CPHF) or Kohn-Sham (CPKS) equations. For uncoupled problems the above density functions can be used in a simple sum-over-states (SoS) framework. For this purpose, instead of the MO-coefficients of only occupied orbitals, also the MO-coefficients of unoccupied orbitals have to be considered, weighted by the orbital energy differences.

III Implementation

The above formalism was implemented within a two-component ZORA versionvan Wüllen (2010) of the quantum chemistry program package TurbomoleAhlrichs et al. (1989). The derivative operator $\hat{\bm{\partial}}(\vec{r}^{\prime}\vee\vec{r})$ was limited to first and second derivatives. The tensor operator can always be reduced to any sum of tensor products of the electronic position vector with respect to some arbitrary origin $\vec{r}-\vec{r}_{\mathrm{O}}$ or with respect to the center of a nucleus $A$: $\vec{r}-\vec{r}_{A}$. Eq. (28) is rewritten in terms of two-component density matrices $\mathbf{D}^{(k)}$:

[TABLE]

Within above mentioned restrictions for the derivative operator, a fortran code was automatically generated for evaluation of the integrand of eq. (29) (term in square brackets) with the computer algebra system Mathematica.Wolfram Research, Inc. (2016)

The code was generated for the derivative operators $\vec{\nabla}\otimes$, $\vec{\nabla}\cdot$, $\vec{\nabla}\times$, $\vec{\nabla}\otimes\vec{\nabla}\cdot$, $\vec{\nabla}\cdot\vec{\nabla}\otimes$, $\vec{\nabla}\otimes\vec{\nabla}\times$, and $\vec{\nabla}\times\vec{\nabla}\times$. For computation of maximally second derivatives, at most third derivatives of basis functions are required due to appearance of first derivatives of matrix elements containing the small component. Furthermore, up to second derivatives of the model potential are needed with respect to $\vec{r}$. These were implemented following Ref. van Wüllen, 1998. In case of second derivatives, third derivatives of density functionals with respect to the model density are required. These were approximated by central finite differences of analytical second derivatives of the density functional.

Within our implementation, integrals are evaluated by default on a grid using standard numerical integration methods, established for density functional theory (DFT) calculations. Integrals that do not contain the ZORA-factor $\omega$ can typically be evaluated analytically. In the present implementation, however, only a few analytical integrals are available.

We added also a SoS module for calculation of second order properties within an uncoupled DFT approach. The SoS calculation is carried out with transition density functions which are evaluated analogously to eqs. (20a) to (20d) as densities between occupied and unoccupied orbitals $i,a$:

[TABLE]

and so on. The second order SoS expression is

[TABLE]

where $\epsilon$ are the orbital energies, $i$ are indices of occupied and $a$ are indices of unoccupied orbitals. The integrals $\Braket{\psi_{i}}{\hat{\mathbf{O}}}{\psi_{a}}$ are evaluated via the corresponding generic transition one-electron tensor functions:

[TABLE]

with

[TABLE]

The contraction of AO-matrix elements with LCAO-coefficients is implemented via density matrices as described in the last section or matrix multiplications of the type $(\mathbf{C}_{\mathrm{occ}}^{(\xi)})^{\dagger}\mathbf{O}^{IJ,j,k,(m,l)}_{\mathrm{AO}}\mathbf{C}_{\mathrm{unocc}}^{(\xi)}$, where $\xi$ can be $\alpha$ or $\beta$, $\mathbf{C}_{\mathrm{occ}}^{(\xi)}$ and $\mathbf{C}_{\mathrm{unocc}}^{(\xi)}$ are the $\alpha$- and $\beta$-spin blocks of the block of occupied and unoccupied orbitals in the coefficient matrix, respectively. Here, for component $(m,l)$ of the matrix of integrals in AO-basis $\mathbf{O}^{IJ,j,k,(m,l)}_{\mathrm{AO}}$ the sum of matrix multiplications with $\alpha$- and $\beta$-coefficient matrices, that correspond to $\bm{\varsigma}(m,l)$ is formed.

In order to account for renormalisation effects, which are in the ZORA approach particularly relevant for the description of the energetically lowest lying orbitals with main contributions close to the nucleus, we implemented the possibility of renormalisation of the ZORA wave function by redefinition of the coefficient matrices:

[TABLE]

This renormalisation can be evaluated directly within the general implementation of one-electron operators presented above.

IV Computational details

Quasi-relativistic two-component calculations are performed within ZORA at the level of complex generalized Hartree–Fock (cGHF) or Kohn–Sham (cGKS) with a modified versionvan Wüllen (2010) of the quantum chemistry program package TurbomoleAhlrichs et al. (1989).

For Kohn–Sham (KS)-DFT calculations 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. In comparison to relativistic coupled cluster (CC) calculations this functional performed well for the description of $\mathcal{P,T}$-odd effects in diatomic radicals in our previous work, which motivates the present choice.Gaul and Berger (2017); Gaul et al. (2019); Gaul and Berger (2018)

For all calculations a 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 by $\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 Tl and Ra. The largest exponent coefficients of the s, p, d and f subsets are $2\times 10^{9}\leavevmode\nobreak\ a_{0}^{-2}$, $5\times 10^{8}\leavevmode\nobreak\ a_{0}^{-2}$, $13300.758\leavevmode\nobreak\ a_{0}^{-2}$ and $751.8368350\leavevmode\nobreak\ a_{0}^{-2}$, respectively. This basis set has proven successful in calculations of nuclear-spin dependent $\mathcal{P}$-violating interactions and $\mathcal{P,T}$-odd effects induced by an eEDM in heavy polar diatomic molecules.Isaev and Berger (2012, 2013, 2014); Gaul and Berger (2017); Gaul et al. (2019); Gaul and Berger (2018) The F atom was represented with a decontracted atomic natural orbital (ANO) basis set of triple-$\zeta$ qualityRoos et al. (2004).

The ZORA-model potential $\tilde{V}(\vec{r})$ as proposed by van Wüllenvan Wüllen (1998) was employed with additional dampingLiu et al. (2002).

For calculations of two-component wave functions and properties a finite nucleus was used, described by a normalized spherical Gaussian nuclear density distribution $\rho_{\mathrm{nuc},A}(\vec{r})=-\frac{\zeta_{A}^{3/2}}{\pi^{3/2}}\mathrm{e}^{-\zeta_{A}\left\lvert\vec{r}-\vec{r}_{A}\right\rvert^{2}}$, where $\zeta_{A}=\frac{3}{2r_{\text{nuc},A}^{2}}$ and the root mean square radius $r_{\text{nuc},A}$ of nucleus $A$ was used as suggested by Visscher and Dyall.Visscher and Dyall (1997) The mass numbers $A$ were chosen to correspond to the isotopes 205Tl,223Ra.

Nuclear equilibrium distances were obtained at the levels of GHF-ZORA and GKS-ZORA/B3LYP, respectively. As convergence criteria an energy change of less than $10^{-5}\leavevmode\nobreak\ E_{\text{h}}$ was used. For DFT calculations of analytic energy gradients with respect to the displacement of the nuclei the nuclei were approximated as point charges. The equilibrium distances obtained are for 205TlF 2.08 Å (cGHF) and 2.12 Å (cGKS) and for 223RaF 2.28 Å (cGHF) and 2.26 Å (cGKS).

All properties were computed with and without renormalisation of the wave function according to (34).

V Results and discussion

V.1 Enhancement of various sources of $\mathcal{P,T}$-violation

in paramagnetic RaF and diamagnetic TlF

A measurement of a permanent electric dipole moment in vanishing electric field would indicate a simultaneous violation of parity $\mathcal{P}$ and time-reversal $\mathcal{T}$ symmetry (see e.g Ref. Khriplovich and Lamoreaux, 1997). Molecular systems provide currently the strictest experimental limits on permanent electric dipole moments.DeMille (2015); Andreev et al. (2018)

There are many possible $\mathcal{P,T}$-violating sources that can lead to a permanent electric dipole moment in a molecule depending on the nuclear and electronic spin states of the molecule.Khriplovich and Lamoreaux (1997); Kozlov and Labzowsky (1995) For paramagnetic molecules such as RaF pronounced sensitivity is expected for a potential permanent electric dipole moment of the electron (eEDM) $d_{\mathrm{e}}$ and from $\mathcal{P,T}$-odd scalar-pseudoscalar nucleon-electron current (SPNEC) interactions $k_{\mathrm{S}}$.Ginges and Flambaum (2004)

However, in case of a non-zero nuclear spin additional contributions from interactions of the electron cloud with a potential permanent electric dipole moment of the proton (pEDM) $d_{\mathrm{p}}$, interactions with a net electric dipole moment of the nucleus, called Schiff moment $\mathcal{S}$, tensor-pseudotensor nucleon-electron current (TPNEC) interactions $k_{\mathrm{T}}$, and pseudoscalar-scalar nucleon-electron current (PSNEC) interactions $k_{\mathrm{p}}$ can occur. The latter, PSNEC, vanishes in the limit of infinitely large mass of the nucleus, which appears to be a reasonable approximation for heavy-elemental molecules, as $m_{\mathrm{e}}\ll m_{\mathrm{nuc}}$.Khriplovich and Lamoreaux (1997)

Nuclear spin-dependent contributions are the dominating sources for $\mathcal{P,T}$-violation in diamagnetic molecules such as TlF, where interactions with an eEDM or due to SPNECs appear only as indirect interactions via hyperfine induced coupling, because of vanishing total effective electron spin.

In paramagnetic molecules containing nuclei with spin quantum number $I$ larger than ${}^{1}/_{2}$, interactions with potential higher $\mathcal{P,T}$-odd nuclear moments, as e.g. nuclear magnetic quadrupole moments (NMQM), could contribute to a permanent molecular EDM.

The full $\mathcal{P,T}$-odd effective spin-rotational Hamiltonian for a paramagnetic diatomic molecule for nucleus $A$ with nuclear spin quantum number larger than ${}^{1}/_{2}$ reads (for parts of ${\hat{H}_{\mathrm{sr}_{A}}}$ see Refs. Kozlov and Labzowsky, 1995 and Hinds and Sandars, 1980)

[TABLE]

Here $\vec{\lambda}$ is the unit vector pointing from the heavy to the light nucleus, $\hat{\vec{S}}$ is the effective electron spin of the molecule, $\hat{\vec{I}}_{A}$ the effective nuclear spin of nucleus $A$, $\hat{\mathbf{T}}_{A}$ is a second-rank tensor which can be constructed from the components of $\hat{\vec{I}}_{A}$ (for details see Ref. Kozlov and Labzowsky, 1995), $\tilde{\mathcal{M}}_{A}=\frac{-1}{2I_{A}(2I_{A}-1)}\mathcal{M}_{A}$ with the NMQM $\mathcal{M}$ (see Ref. Kozlov and Labzowsky, 1995), $\Omega$ is the projection of the effective electron spin on the molecular axis and $\mathcal{I}_{A}$ is the projection of $\hat{\vec{I}}_{A}$ on the molecular axis. For a diamagnetic molecule all terms that depend on the effective electronic spin $\vec{S}$ (here all proportional to $\Omega$ and $\Theta$) vanish. The constants $W$ are electronic structure coupling constants enhancing $\mathcal{P,T}$-violating parameters in molecules, that need to be determined by electronic structure calculations. $R_{\text{vol}}$ is a nuclear structure factor that enhances the pEDM and can be determined from nuclear structure calculations.

In the following we will focus on the electronic structure enhancement factors $W$ in RaF and TlF. Thereby we will not calculate $W_{\mathrm{p}_{A}}$ which is supposed to be many orders of magnitude smaller than the other effects.Khriplovich and Lamoreaux (1997); Ginges and Flambaum (2004) Furthermore, we do not include in our discussion those nuclear-spin dependent effects $W^{\mathrm{m}}_{\mathrm{s}_{A}}$ and $W^{\mathrm{m}}_{\mathrm{d}_{A}}$ that are induced by hyperfine coupling and thus are second-order molecular properties that can be obtained from a linear response treatment. Assuming the molecular axes to be aligned along the $z$-axis the remaining electronic structure parameters are defined in the following (see e.g Khriplovich and Lamoreaux, 1997; Kozlov and Labzowsky, 1995; Mårtensson-Pendrill and Öster, 1987; Hinds and Sandars, 1980; Quiney et al., 1998):

i) electronic structure enhancement of the eEDM

[TABLE]

with $e$ being the elementary charge and $\hbar=\frac{h}{2\pi}$ being the reduced Planck’s constants;

ii) electronic structure enhancement of the NMQM

[TABLE]

with the constant $k_{\mathrm{em}}$ being $\frac{\mu_{0}}{4\pi}$ in SI units with $\mu_{0}$ being the magnetic constant (see Ref. Cohen et al., 2008 for other choices of $k_{\mathrm{em}}$ that correspond to different unit systems);

iii) electronic structure enhancement of nuclear Schiff moment and volume effect due to a pEDM

[TABLE]

in calculations of $W_{\mathcal{S},A}$ the value was calculated at eight points at a distance of $1.7\times 10^{-25}\leavevmode\nobreak\ a_{0}$ around the nucleus (cubic arrangement) and averaged. This Hamiltonian is a consequence of Schiff’s theoremSchiff (1963) and is the dipole contribution of an expansion of the electric potential of a finite nucleus (for details see Ref. Sushkov, Flambaum, and Khriplovich, 1984);

iv) electronic structure enhancement of the pEDM due to magnetic fields of moving electrons

[TABLE]

with constant $k$, which is 1 in SI units and $c^{-1}$ in Gauss units, the orbital angular momentum operator with respect to nucleus $A$ $\hat{\vec{\ell}}_{A}=(\vec{r}-\vec{r}_{A})\times\hat{\vec{p}}$, the nuclear magneton $\mu_{\mathrm{N}}=\frac{e\hbar}{2m_{\mathrm{p}}}$ with the mass of the proton $m_{\mathrm{p}}$, the nuclear magnetic moment $\mu_{A}$, the nuclear charge $Z_{A}$ and mass number $A_{A}$ of nucleus $A$. Here we assume that the contributions stem from a single active valence proton in the nuclear shell (see Ref. Hinds and Sandars, 1980);

v) electronic structure enhancement of SPNEC interactions

[TABLE]

with the Fermi weak coupling constant $G_{\mathrm{F}}=2.22249\times 10^{-14}E_{\mathrm{h}}a_{0}^{3}$;

vi) electronic structure enhancement of TPNEC interactions

[TABLE]

With the method presented in the previous section all these properties can be directly evaluated as illustrated in Table 1 on page 1. We choose as exemplary molecular systems RaF and TlF, which were well studied before and are considered promising candidates for experiments that aim at a measurement of $\mathcal{P,T}$-odd properties.Isaev and Berger (2013); Kudashov et al. (2014); Sasmal et al. (2016); Hinds and Sandars (1980); Parpia (1997); Quiney et al. (1998)

We neglect any magnetic or many-electron effects on the above presented properties that may arise from the ZORA transformation as these are expected to be low for heavy elements. For some of the properties such effects were discussed elsewhere (see e.g. Refs. Lindroth, Lynn, and Sandars, 1989; Gaul et al., 2019).

V.1.1 Diamagnetic Molecules: TlF

In this section we discuss results for TlF that were obtained from property calculations with an implementation of the approach described in the previous sections. Results of our calculations on the level of cGHF- and cGKS-B3LYP-ZORA are compared to results from literature in Table 2 on page 2.

We see an excellent agreement (deviations $\leq 5\leavevmode\nobreak\ \%$) between cGHF-ZORA and DHF results reported by Quiney et. al. for all calculated properties.Quiney et al. (1998) Furthermore, renormalisation of the wave function does not play an important role for the nuclear spin-dependent properties and is always below 1 %. This is to be expected as all major contributions stem from the valence molecular orbitals at the position of the nucleus. This shows that ZORA is appropriate for the quantitative description of relativistic effects due to nuclear-spin dependent $\mathcal{P,T}$-odd interactions in a heavy molecule, which reinforces related previous findings for nuclear-spin independent and nuclear-spin dependent $\mathcal{P}$-odd interactions in molecules Berger, Langermann, and van Wüllen (2005); Berger and van Wüllen (2005); Berger and Stuber (2007); Nahrwold and Berger (2009); Isaev, Hoekstra, and Berger (2010); Isaev and Berger (2012, 2014). Correlation effects were estimated on the DFT level with the B3LYP hybrid density functional. This functional, however, seems to overshoot electron correlation effects leading to values being too low in magnitude with deviations of up to $23\leavevmode\nobreak\ \%$ in comparison to GRECP/RCC-SD calculations.Petrov et al. (2002)

V.1.2 Paramagnetic Molecules: RaF

We report results on electronic structure enhancement factors of nuclear spin-dependent and nuclear spin-independent $\mathcal{P,T}$-effects in 223RaF in Table 3 on page 3 and compare to all available literature data.

Also for nuclear-spin independent properties and the nuclear- and electron-spin dependent NMQM, most important contributions stem from valence molecular orbitals and thus effects of renormalisation are negligible. The reasonable agreement between cGHF/cGKS-ZORA and four component coupled cluster calculationsKudashov et al. (2014); Sasmal et al. (2016) for $W_{\mathrm{d}}$ and $W_{\mathrm{s}}$ was discussed elsewhere.Kudashov et al. (2014); Gaul and Berger (2017); Gaul et al. (2019) Our present values for $W_{\mathrm{d}}$ and $W_{\mathrm{s}}$ differ in the last reported digit from the results in Refs. Gaul and Berger, 2017; Gaul et al., 2019 as we consider herein a different isotope of Ra. Our calculations of $W_{\mathcal{S}}$ are with deviations of $8\leavevmode\nobreak\ \%$ (cGHF) in good agreement with GRECP-FSCC calculations (Ref. Kudashov et al., 2014). As observed in previous studies,Gaul and Berger (2017); Gaul et al. (2019) DFT tends to give results too low in absolute value. However, the value of $W_{\mathcal{S}}$ still agrees reasonably (deviation is $12\leavevmode\nobreak\ \%$) with the coupled cluster calculations. From this we expect a similar precision of cGHF- and cGKS-ZORA calculations of $W_{\mathrm{T}}$ and $W_{\mathrm{m}}$.

Our calculations show that nuclear-spin dependent $\mathcal{P,T}$-odd enhancement factors in RaF are by about a factor of $1/2$ to $1/8$ smaller in magnitude compared to TlF. The ratios of $\frac{W_{\mathrm{T}}}{W_{\mathcal{S}}}$ and $\frac{W_{\mathrm{m}}}{W_{\mathcal{S}}}$ are considerably different for TlF and RaF. Thus data from measurements of RaF would complement data from TlF measurements as different regions in the parameter space of the $\mathcal{P,T}$-odd parameters $k_{\mathrm{T}}$, $d_{\mathrm{p}}$ and $\mathcal{S}$ are covered.

In Table 3 on page 3 we presented values for $W_{\mathcal{M}}$, as well. In comparison to predictions of NMQM enhancement made for other moleculesFlambaum, DeMille, and Kozlov (2014) the values for RaF are very large (as large as predicted for YbF or ThO). Furthermore the 223Ra nucleus is known to have a octupole deformation, which is expected to enhance NMQM effects significantly on the nuclear structure level.Flambaum (1994) This makes 223RaF a promising candidate for setting strict limits on NMQM induced permanent electric dipole moments.

VI Conclusion

We outlined a generally applicable approach to the evaluation of arbitrary relativistic properties within an approximate quasi-relativistic wave function. Automated code generation via a computer algebra system was applied to obtain a pilot implementation within a modified version of a quantum chemical program for ZORA wave functions. Within the approach presented herein, relativistic first and second order properties ranging from commonly available molecular properties, such as NMR-shielding constants, to less common discrete symmetry violating properties are accessible in a single implementation. The flexibility of this property toolbox was demonstrated by computation of a number of $\mathcal{P,T}$-odd effects that are important for fundamental physics research with diatomic molecules

Within this study large enhancements of the nuclear magnetic quadrupole moment and nuclear spin-dependent as well as nuclear spin-independent sources of $\mathcal{P,T}$-violation in RaF were determined. This shows that 223RaF is a well suited system for setting strict limits on $\mathcal{CP}$-violation in essentially all sectors of particle physics.

Acknowledgements.

Computer time provided by the center for scientific computing (CSC) Frankfurt and financial support by the Deutsche Forschungsgemeinschaft via Sonderforschungsbereich 1319 (ELCH) “Extreme Light for Sensing and Driving Molecular Chirality” are gratefully acknowledged.