Significance of Non-Linear Terms in the Relativistic Coupled-Cluster Theory in the Determination of Molecular Properties

Abstract

The relativistic coupled-cluster (RCC) theory has been applied recently to a number of heavy molecules to determine their properties very accurately. Since it demands large computational resources, the method is often approximated to single and double excitations (RCCSD method). The effective electric fields (${\mathcal{E}}_{eff}$) and molecular permanent electric dipole moments (PDMs) of SrF, BaF, and mercury monohalides (HgX with X = F, Cl, Br, and I) molecules are of immense interest for probing fundamental physics. In our earlier calculations of ${\mathcal{E}}_{eff}$ and PDMs for the above molecules, we neglected the non-linear terms in the property evaluation expression of the RCCSD method. In this work, we demonstrate the roles of these terms in determining the above quantities and their computational time scalability with the number of processors of a computer. We also compare our results with previous calculations that employed variants of RCC theory, as well as other many-body methods and available experimental values.

1. Introduction

The coupled-cluster (CC) theory is considered to be the gold standard of electronic structure calculations in atoms and molecules [1,2]. It owes the title to its ability to capture electron correlation effects to a much better extent than other well-known many-body approaches such as configuration interaction (CI) [3], at a given level of truncation. This feature has led to accurate calculations of many properties in both the atomic and molecular systems (for example, see [4,5]). We shall focus on the application of this method to evaluate molecular properties that are useful to probe fundamental physics, specifically the permanent electric dipole moment (PDM) and parity and time-reversal violating electric dipole moment of the electron (eEDM) [6,7]. The molecular PDM is a very interesting property, and it plays a role in the sensitivity of an eEDM experiment through the polarizing factor [8,9]. The PDM is also an extremely relevant property in the ultracold sector, and molecules with large PDMs find innumerable applications in that domain. For example, the SrF molecule possesses a fairly large PDM and hence gives rise to long-range, tunable, and anisotropic dipole-dipole interactions. This aspect, in combination with the fact that SrF is laser-coolable, makes the molecule important for applications such as exploring new quantum phases and quantum computing [10].

The extremely tiny eEDM is yet to be detected. The upper bounds to it are extracted by a combination of relativistic many-body theory and experiment [11]. These bounds, in turn, help to constrain several theories that lie beyond the Standard Model of particle physics, for example supersymmetric theories [12]. The knowledge of the eEDM also aids in understanding the underlying physics that describes the matter-antimatter asymmetry in the universe [13]. The theoretical molecular property of interest to eEDM is the effective electric field, ${\mathcal{E}}_{\mathrm{eff}}$. It is the internal electric field that is experienced by an electron due to other electrons and nuclei in a molecule. An accurate estimate of this quantity is used in setting or improving an upper bound to eEDM (for example, [5]), or to propose a new candidate for molecular eEDM experiments (for example, [8]). This quantity can only be obtained using a relativistic many-body theory [14]. Calculating the PDM provides information on the polarizing factor for molecules that are proposed for eEDM searches, where the property has not been measured.

There have been several calculations of ${\mathcal{E}}_{\mathrm{eff}}$ for various molecules using the single and double excitations approximation in the relativistic CC theory (RCCSD method), for example [15,16]. In our earlier RCCSD calculations [5,8,17,18,19,20], the expectation value evaluating the expression was approximated to only the linear terms (referred to as LERCCSDmethod). Later, the calculations performed for HgX (X = F, Cl, Br, and I), SrF, and BaF besides other molecules were verified by using the finite-field energy derivative approach of the RCCSD theory (FFRCCSD method) [9], by adding the interaction Hamiltonians along with the residual Coulomb interaction operator. The LERCCSD and the FFRCCSD approaches showed excellent agreements (within one percent) in the values of ${\mathcal{E}}_{\mathrm{eff}}$. The results for the PDMs obtained in these methods were comparable for SrF and BaF and also overestimated the property with respect to their experimental values, but they differed substantially for HgX (with as much as 20 percent for HgI). The shortcomings of the above FFRCCSD method were that the accuracy of the results depended on numerical differentiation. Moreover, orbital relaxation effects were neglected by not including the perturbation in the Dirac–Fock (DF) level itself, in order to avoid breaking of Kramer’s symmetry in the presence of a time-reversal symmetry violating eEDM interaction, which has to be compensated for eventually with further iterations.

Here, we intend to calculate the values of ${\mathcal{E}}_{eff}$ and PDM by including the non-linear terms in the expectation value evaluation expression of the RCCSD method (nLERCCSD method). We adopt the intermediate-diagram approach as discussed in [21,22] to implement these non-linear RCC terms. For this purpose, we undertook molecules that were very relevant for eEDM studies. HgX molecules were identified as promising candidates for future eEDM searches, owing to their extremely large effective electric fields, as well as experimental advantages [8]. A recent work that proposed to laser-cool HgF opened new avenues for an upcoming eEDM experiment with the molecule [23]. Another very important molecule in this regard is BaF, and two eEDM experiments are simultaneously underway for this system [24,25]. Experimental values of the PDMs are available only for BaF among the systems that we mentioned above. We also present results for the PDM of SrF, as it was the first molecule to be laser-cooled [10], and a very precise measurement of this quantity has been reported [26].

2. Theory and Implementation

In the RCC theory, the wave function of a molecular state is expressed as [27]:

$$\begin{array}{c}\hfill |\mathsf{\Psi }\rangle =e^T|{\mathsf{\Phi }}_0\rangle ,\end{array}$$

(1)

where T is the cluster operator and $|{\mathsf{\Phi }}_0\rangle$ is the reference state obtained by mean-field theory. We used the Dirac–Coulomb Hamiltonian in our calculations, and $|{\mathsf{\Phi }}_0\rangle$ was obtained using the DF method. In the RCCSD method, we approximated $T=T_1+T_2$ with subscripts 1 and 2 indicating single and double excitations, respectively, and they are given using the second-quantization operators as:

$$\begin{array}{ccc}\hfill T_1& =& \sum _{i,a}{t}_i^a{a}_a^{†}{a}_i\mathrm{and}{T}_2=\frac{1}{2}\sum _{i,j,a,b}{t}_{ij}^{ab}{a}_a^{†}{a}_b^{†}{a}_j{a}_i,\hfill \end{array}$$

(2)

where the notation $i,j$ is used to denote holes, $a,b$ refer to particles, $t_i^a$ is the one hole-one particle excitation amplitude, and $t_{ij}^{ab}$ is the two-hole two-particle excitation amplitude.

We employed the UTChem [28,29,30] program for the DF calculations, the atomic orbital to molecular orbital integral transformations, as well as generating the property integrals and the Dirac08 [31] code to obtain the RCCSD excitation operator amplitudes. It is important to reiterate that all the non-linear terms were included in the equations of the RCCSD method to determine the excitation amplitudes.

The expectation value of an operator, O, in the (R)CC method, can be written as follows:

$$\begin{array}{ccc}\hfill \langle O\rangle =\frac{\langle \mathsf{\Psi }\left|O\right|\mathsf{\Psi }\rangle }{\langle \mathsf{\Psi }|\mathsf{\Psi }\rangle }& =& \frac{\langle {\mathsf{\Phi }}_0\left|e^{T^{†}}O{e}^T\right|{\mathsf{\Phi }}_0\rangle }{\langle {\mathsf{\Phi }}_0\left|e^{T^{†}}{e}^T\right|{\mathsf{\Phi }}_0\rangle }\hfill \\ & =& {\langle {\mathsf{\Phi }}_0\left|\left(e^{T†}O{e}^T\right)\right|{\mathsf{\Phi }}_0\rangle }_c,\hfill \end{array}$$

(3)

where the subscript, “c”, means that each term is fully contracted [32], or in the diagrammatic terminology, connected [22].

The PDM of a molecule is determined as [33]:

$$\begin{array}{c}\hfill \mu =\langle \mathsf{\Psi }\left|D\right|\mathsf{\Psi }\rangle +\sum _A{Z}_A{r}_A,\end{array}$$

(4)

where D is the electric dipole operator, the index A runs over the number of nuclei, $Z_A$ is the atomic number of the Ath nucleus, and $r_A$ is the position vector from the origin to the site of the Ath nucleus. The first term in the above expression is the electronic term, while the second term is the nuclear contribution.

Similarly, ${\mathcal{E}}_{\mathrm{eff}}$ is evaluated as:

$$\begin{array}{ccc}\hfill {\mathcal{E}}_{\mathrm{eff}}& =& \langle \mathsf{\Psi }|\sum _{i=1}^{N_e}\beta {\Sigma }_i·E_i^{\mathrm{intl}}|\mathsf{\Psi }\rangle ,\hfill \end{array}$$

(5)

where the summation is over the number of electrons, $N_e$, $\beta$ is one of the Dirac matrices (also known as ${\gamma }_0$ in the literature), $\Sigma$ is the (4 × 4) version of Pauli matrices, and $E_i^{\mathrm{intl}}$ is the internal electric field that is experienced by the ith electron, and is given by the negative of the gradient of the sum of electron-nucleus and electron-electron interaction potentials. Since the expression given above involves evaluating integrals over a two-body Coulomb operator, $\frac{1}{r_{ij}}$, and is complicated, we resorted to employing an effective eEDM Hamiltonian instead of the one introduced above [14]. It follows that:

$$\begin{array}{c}\hfill {\mathcal{E}}_{\mathrm{eff}}=-2ic\langle \mathsf{\Psi }|\sum _{i=1}^{N_e}\beta {\gamma }_5{p}_i^2|\mathsf{\Psi }\rangle ,\end{array}$$

(6)

where ${\gamma }_5$ is the product of the gamma matrices (given by $i{\gamma }_0{\gamma }_1{\gamma }_2{\gamma }_3$), while $p_i$ is the momentum of the ith electron.

In the LERCCSD method, the following expression has been used in the evaluation of the expectation values:

$$\begin{array}{c}\hfill \langle O\rangle ={\langle {\mathsf{\Phi }}_0\left|{(1+T_1+T_2)}^{†}O(1+T_1+T_2)\right|{\mathsf{\Phi }}_0\rangle }_c.\end{array}$$

(7)

These terms are represented using Goldstone diagrams and are shown in Figure 1. Note that $O{T}_2$ and its Hermitian conjugate are zero, due to Slater–Condon rules [34,35]. Diagrammatically, such a diagram will have at least two open lines, that is it is not fully connected. The evaluation of the properties using the LERCCSD approximation misses the contributions corresponding to many correlation effects that will arise from the relativistic third-order many-body perturbation theory (RMBPT method). On the other hand, it is not possible to evaluate exactly Equation (3) even in the RCCSD method as it contains a non-terminating expression. However, it is possible to demonstrate the importance of contributions from the leading-order non-linear terms corresponding to the third- and fourth-order effects of the RMBPT method. It is still extremely challenging to perform direct calculations by incorporating the non-linear terms of Equation (3) in heavier molecules, due to the amount of computations involved in it. In order to tackle this issue, we adopted an additional computational step by breaking the non-linear terms into intermediate parts as described more elaborately below. Further, we parallelized the program using Message Passing Interface (MPI) and show the scalability of their calculations with the number of processors of a computer.

The approach can be understood by revisiting the diagrams in Figure 1, and relating each of them to Figure 2, Figure 3, and Figure 4. As an example, we consider Sub-figure (ii). The property vertex has one incoming particle line and an outgoing hole line. We define it as a particle-hole vertex. Such particle-hole vertices can be found in Sub-figures (v) to (viii) as well. In the intermediate-diagram formalism, the vertex O is removed and replaced successively by each of the particle-hole (p-h) diagrams (more precisely, their Hermitian conjugates) of Figure 4. We assign the notation $O_{p-h}$ for such a vertex. This sequence of operations already generates 28 diagrams from $O_{p-h}{T}_1$ and includes terms that occur in the RMBPT method. We note that $O_{p-h}{T}_1$ in the case of the Hermitian conjugate of Sub-figure (i) of Figure 4 gives back the LERCCSD diagram for $O{T}_1$. Similar to $O_{p-h}$, we also construct the analogous $O_{h-h}$ and $O_{p-p}$ diagrams (as given in Figure 2 and Figure 3, respectively), and generate more terms. We note that the property vertex from the DF diagram in Figure 1 is not replaced with any intermediate diagram, as otherwise there would be repetition of diagrams in the calculations. Further, one has to be careful to avoid repetition of diagrams while contracting effective O operators with the T RCC operators. For example, it can be shown that $T_1^{†}O{T}_1$ diagrams can appear twice through $T_1^{†}{O}_{p-h}$ and $T_1^{†}{O}_{p-p/h-h}{T}_1$ terms. Such diagrams are identified by careful analysis, and their double counting is removed manually.

As can be seen from the above discussions, some of the diagrams that were undertaken in this procedure demanded up to the order of $n_h^3{n}_p^3$ in computational cost for $n_h$ number of holes and $n_p$ number of particles. Therefore, the intermediate diagram approach systematically took into account non-linear terms while simultaneously cutting down drastically on the computational cost as compared to a direct brute-force evaluation of a non-linear expectation value expression. This could be understood by choosing an example as follows. Replacing O of Figure 1v by the property vertex with Sub-figure (xxv) from Figure 4 entailed a computational cost $\mathcal{O}\left(n_p^4{n}_h^4\right)$ for the direct evaluation of such a diagram. However, the intermediate-diagram approach led to a cost of $\mathcal{O}(n_p^2{n}_h^2+n_p^3{n}_h^3)$. This became especially relevant when we performed computations on heavy systems and with high-quality basis sets, such as those that we chose for this work. For instance, the RCC calculations on HgF involved $n_h=89$ and $n_p=429$, and therefore, the computational cost with an intermediate-diagram approach was a full five orders of magnitude smaller than a brute-force approach to computing the same diagram (without considering any molecular point group symmetries). A similar level of reduction in computational cost could be seen from the heaviest HgI as well. We add at this point that we exploited the $C_8$ double point group symmetry in our nLERCCSD code, as we had done for the earlier LERCCSD program [5,36]. This aspect also substantially lessened the computational cost, as it restricted the number of matrix elements to be evaluated based on group theoretic considerations. For example, $O{T}_1$ involves computing matrix elements of the form $\langle a\left|O\right|i\rangle$. Given that we had 89 holes and 429 particles, the number of possible matrix elements was ∼$3.8\times {10}^5$. However, since we imposed the restriction that both i and a should belong to the same irreducible representation, we needed to evaluate only ∼$7.2\times {10}^4$ matrix elements. Similar considerations for the more complicated terms involving $T_2$ led to evaluating much fewer matrix elements.

3. Results and Discussions

To carry out the calculations in the considered molecules, we chose values for the bond lengths as 2.00686 $A^{\circ }$, 2.42 $A^{\circ }$, 2.62 $A^{\circ }$, 2.81 $A^{\circ }$, 2.075 $A^{\circ }$, and 2.16 $A^{\circ }$ for HgF, HgCl, HgBr, HgI, SrF, and BaF, respectively [37,38,39,40,41,42]. It is to be noted that the chosen values for the HgX molecules were from theory, while those for SrF and BaF were from experiment. Furthermore, we opted for Dyall’s quadruple zeta (QZ) basis for Hg and I [43], Dunning’s correlation consistent polarized valence quadruple zeta (cc-pVQZ) basis for the halide atoms (F, Cl, and Br) [44], and Dyall’s QZ functions augmented with Sapporo’s diffuse functions [45] for Sr and Ba. We chose Dyall’s basis for Hg and I as it was among the most reliable and widely used basis functions for heavy atoms. We did not add diffuse functions as it increased the computational cost drastically for QZ quality basis sets. Moreover, it was found that the inclusion of diffuse functions changed the effective electric field by around 2.5 percent for HgF, and it was expected to lead to a similar difference for the heavier HgX systems [9]. However, in the foreseeable future, such computations could be performed to improve the calculated values of the PDMs. To minimize steep computational costs that we incurred due to our choice of QZ basis sets, as well as performing all-electron calculations, we cut-off the high-lying virtuals above 1000 atomic units (a.u.) for HgX and BaF. At such a high cut-off value, we could expect that the missing contributions would be minimal, and possibly even negligible.

In

Table 1. Contributions from the Dirac–Fock (DF), LERCCSD, and nLERCCSDmethods to the ${\mathcal{E}}_{\mathrm{eff}}$s (in GV/cm) and permanent electric dipole moments (PDMs) (in Debye) of HgX, SrF, and BaF molecules from the present work (denoted as “This work” in the table). Comparison of the two properties from various works with our results are also presented. CASSCF, complete active space self-consistent field; RSPT2, second-order Rayleigh–Schrodinger perturbation theory; X2C-FSCC, exact two-component Hamiltonian–Fock space coupled-cluster.

Molecule	Method	PDM	${\mathcal{E}}_{\mathbf{eff}}$
SrF	CASSCF-MRCI [46]	3.36
	CASSCF-RSPT2 [46]	3.61
	Z-vector [47]	3.45
	LERCCSD [9,33]	3.6	2.17
	FFCCSD [9]	3.62	2.16
	X2C-MRCI [48]	3.20
	X2C-FSCC [48]	3.46
	DF (This work)	2.99	1.54
	LERCCSD (This work )	3.57	2.15
	nLERCCSD (This work)	3.60	2.16
	Experiment [26]	3.4676(1)
BaF	MRCI [49]	2.96
	LERCCSD [33]	3.4	6.50
	FFCCSD [9]	3.41	6.46
	X2C-MRCI [48]	2.90
	X2C-FSCC [48]	3.23
	Z-vector [50]	3.08
	ECP-RASSCF [51]		7.5
	RASCI [52]		7.28
	MRCI [53]		5.1
	MRCI [54]		6.1
	DF (This work)	2.61	4.81
	LERCCSD (This work)	3.32	6.45
	nLERCCSD (This work)	3.37	6.39
	Experiment (PDM) [55]	3.17(3)
HgF	CI [56]	4.15	99.26
	LERCCSD [57]	2.61
	MRCI [53]		68
	MRCI [54]		95
	DF (This work)	4.11	105.69
	LERCCSD [8]		115.42
	FFCCSD [9]	2.92	116.37
	LERCCSD (This work)	3.25	114.93
	nLERCCSD (This work)	3.45	113.77
HgCl	CI [58]	3.28
	LERCCSD [57]	2.72
	LERCCSD [8]		113.56
	FFCCSD [9]	2.96	114.31
	DF (This work)	4.30	104.33
	LERCCSD (This work)	3.26	112.51
	nLERCCSD (This work)	3.45	110.94
HgBr	CI [58]	2.62
	LERCCSD [57]	2.36
	LERCCSD [8]		109.29
	FFCCSD [9]	2.71	109.56
	DF (This work)	4.14	99.72
	LERCCSD (This work)	2.62	109.38
	nLERCCSD (This work)	2.94	107.42
HgI	LERCCSD [57]	1.64
	LERCCSD [8]		109.3
	FFCCSD [9]	2.06	109.56
	DF (This work)	3.61	99.27
	LERCCSD (This work)	1.50	110.00
	nLERCCSD (This work)	2.01	107.38

* The bond lengths chosen in our work were 2.00686 $A^{\circ }$, 2.42 $A^{\circ }$, 2.62 $A^{\circ }$, 2.81 $A^{\circ }$, 2.075 $A^{\circ }$, and 2.16 $A^{\circ }$ for HgF, HgCl, HgBr, HgI, SrF, and BaF, respectively. We used Dyall’s quadruple zeta (QZ) basis for Hg and I, Dunning’s correlation consistent polarized valence quadruple zeta (cc-pVQZ) basis for the halide atoms (F, Cl, and Br), and Dyall’s QZ functions augmented with Sapporo’s diffuse functions for Sr and Ba.

, we present our results for HgX, SrF, and BaF, all using QZ basis sets. We discuss the trends in the PDMs and ${\mathcal{E}}_{\mathrm{eff}}$s across HgX in the three different approaches, namely the LERCCSD and nLERCCSD methods, while briefly making a comparison with the FFRCCSD method from [9] wherever relevant, and we also examined the correlation effects in the property from lower to all-order methods. SrF and BaF molecules were treated as stand-alone systems. Firstly, we observed that the effect of non-linear corrections was to increase the PDM and decrease the effective electric field (except in the case of SrF, where the difference was still within 0.5 percent). We find from Table 1 that for SrF and BaF, the nLERCCSD method yielded PDMs that were very close to their LERCCSD counterparts (within 1.5 percent of each other for both the molecules), but were not in better agreement with experiments than their LERCCSD counterparts. However, the nLERCCSD values agreed well with the results from the earlier work that used the FFRCCSD approach (within 1.2 percent of each other) that also employed a QZ quality basis with diffuse functions. Such a comparison could not be made with the HgX molecules, as available FFRCCSD data used a double zeta (DZ) quality basis. For HgX systems, we observed that unlike in the cases of SrF and BaF, the difference between the LERCCSD and the nLERCCSD results widened from about six percent for HgF and HgCl, to about 25 percent for HgI. The values for ${\mathcal{E}}_{\mathrm{eff}}$ for SrF and BaF showed that the LERCCSD, nLERCCSD, and FFRCCSD methods all agreed to within one percent. In the case of HgX molecules, the LERCCSD and nLERCCSD results were found to differ by at most 2.5 percent. We chose HgF as a representative molecule and performed FFRCCSD calculations with a QZ basis, and we found that its effective electric field was 110.87 GV/cm, which was lesser than the nLERCCSD value by 2.5 percent.

The individual contributions that arose from the diagrams given in Figure 1 to the effective electric fields and PDMs of HgX, SrF, and BaF molecules are given in

Table 2. Individual correlation contributions to the effective electric fields (in GV/cm) of mercury monohalides (HgX; X = F, Cl, Br, and I), SrF, and BaF, from the LERCCSD (abbreviated as “L”) and nLERCCSD (denoted by “nL”) methods. In the first column, A could be O (which corresponds to LERCCSD diagrams) or $O_{x-y}$ (which is associated with nLERCCSD diagrams), where “x” and “y” could stand for the corresponding particle or hole line for a given term. The values are all rounded-off to two decimal places for HgX, while numbers that are extremely small in the case of SrF and BaF are denoted in the scientific notation instead.

Molecule		HgF		HgCl		HgBr		HgI		SrF		BaF
Term	Diagram	L	nL	L	nL	L	nL	L	nL	L	nL	L	nL
DF	Figure 1i	105.69		104.33		99.72		99.27		1.54		4.81
$A{T}_1$	Figure 1ii	17.09	13.11	17.05	12.21	19.83	14.76	23.85	15.71	0.63	0.61	1.79	1.60
$T_1^{†}A{T}_1$	Figure 1iii	−1.85	−0.28	−2.01	−0.25	−2.65	−0.62	−3.66	−0.41	−1.86$\times {10}^{-2}$	−1.00$\times {10}^{-3}$	−7.65$\times {10}^{-2}$	−2.00$\times {10}^{-4}$
	Figure 1iv	−1.41	0.16	−1.40	0.28	−1.21	0.47	−1.56	1.16	−9.01$\times {10}^{-3}$	4.80$\times {10}^{-4}$	−6.47$\times {10}^{-2}$	7.60$\times {10}^{-3}$
$T_1^{†}A{T}_2$	Figure 1v	1.19	0.93	0.65	0.29	0.38	−0.11	0.38	−0.27	2.73$\times {10}^{-3}$	1.02$\times {10}^{-3}$	9.46$\times {10}^{-3}$	2.51$\times {10}^{-3}$
	Figure 1vi	0.05	0.08	0.06	0.05	−0.01	−0.07	−0.03	−0.09	−4.91$\times {10}^{-4}$	−7.93$\times {10}^{-4}$	−1.49$\times {10}^{-3}$	−2.39$\times {10}^{-3}$
	Figure 1vii	0.61	0.58	0.92	0.85	0.66	0.32	0.57	0.19	1.43$\times {10}^{-2}$	1.48$\times {10}^{-2}$	7.04$\times {10}^{-2}$	7.13$\times {10}^{-2}$
	Figure 1viii	−1.31	−1.27	−1.24	−1.18	−0.91	−0.63	−1.26	−0.98	9.63$\times {10}^{-3}$	9.91$\times {10}^{-3}$	−2.49$\times {10}^{-2}$	−2.32$\times {10}^{-2}$
$T_2^{†}A{T}_2$	Figure 1ix	−2.50	−2.46	−2.54	−2.49	−2.68	−2.65	−2.93	−2.89	8.58$\times {10}^{-3}$	6.15$\times {10}^{-3}$	3.22$\times {10}^{-2}$	2.17$\times {10}^{-2}$
	Figure 1x	−0.17	−0.17	−0.15	−0.14	−0.14	−0.13	−0.13	−0.11	−2.17$\times {10}^{-3}$	−1.93$\times {10}^{-3}$	−6.87$\times {10}^{-3}$	−6.83$\times {10}^{-3}$
	Figure 1xi	−1.22	−1.40	−1.50	−1.47	−1.65	−1.85	−1.96	−1.99	−1.96$\times {10}^{-2}$	−2.17$\times {10}^{-2}$	−7.54$\times {10}^{-2}$	−7.71$\times {10}^{-2}$
	Figure 1xii	−0.17	−0.17	−0.15	−0.14	−0.14	−0.13	−0.13	−0.11	−2.17$\times {10}^{-3}$	−1.93$\times {10}^{-3}$	−6.87$\times {10}^{-3}$	−6.83$\times {10}^{-3}$
	Figure 1xiii	−1.64	−1.57	−1.67	−1.57	−1.70	−1.58	−1.84	−1.69	−1.38$\times {10}^{-3}$	−1.96$\times {10}^{-3}$	−3.20$\times {10}^{-4}$	−9.61$\times {10}^{-4}$
	Figure 1xiv	−0.10	−0.10	−0.10	−0.10	−0.10	−0.10	−0.10	−0.10	−5.39$\times {10}^{-4}$	−5.53$\times {10}^{-4}$	−2.28$\times {10}^{-3}$	−2.33$\times {10}^{-3}$
	Figure 1xv	0.77	0.74	0.36	0.37	0.08	0.12	−0.37	−0.21	4.42$\times {10}^{-3}$	4.51$\times {10}^{-3}$	2.82$\times {10}^{-3}$	3.18$\times {10}^{-3}$
	Figure 1xvi	−0.10	−0.10	−0.10	−0.10	−0.10	−0.10	−0.10	−0.10	−5.39$\times {10}^{-4}$	−5.53$\times {10}^{-4}$	−2.28$\times {10}^{-3}$	−2.33$\times {10}^{-3}$
Total		114.93	113.77	112.51	110.94	109.38	107.42	110.00	107.38	2.15	2.16	6.45	6.39

and

Table 3. Correlation contributions to the PDMs (in Debye) of mercury monohalides (HgX; X = F, Cl, Br, and I), SrF, and BaF. The notation is the same as in Table 2. The entry “NC” stands for nuclear contribution to the PDM.

Molecule		HgF		HgCl		HgBr		HgI		SrF		BaF
Term	Diagram	L	nL	L	nL	L	nL	L	nL	L	nL	L	nL
DF	Figure 1i	−767.04		−925.61		−1002.31		−1075.83		−375.75		−578.39
$A{T}_1$	Figure 1ii	−0.60	−0.78	−0.83	−1.01	−1.26	−1.54	−1.92	−2.33	0.63	0.65	0.80	0.82
$T_1^{†}A{T}_1$	Figure 1iii	0.21	0.04	0.26	0.06	0.34	0.11	0.48	0.23	0.14	−0.01	0.19	−0.02
	Figure 1iv	−0.45	0.05	−0.48	0.07	−0.62	0.13	−0.79	0.26	−0.18	−0.01	−0.23	−0.02
$T_1^{†}A{T}_2$	Figure 1v	0.10	0.11	0.13	0.13	0.20	0.19	0.30	0.29	2.44$\times {10}^{-2}$	2.53$\times {10}^{-2}$	3.04$\times {10}^{-2}$	3.13$\times {10}^{-2}$
	Figure 1vi	0.01	0.01	0.00	0.00	0.00	0.00	0.00	0.00	−1.99$\times {10}^{-3}$	2.13$\times {10}^{-3}$	2.40$\times {10}^{-3}$	2.55$\times {10}^{-3}$
	Figure 1vii	0.01	0.01	0.01	−0.01	−0.01	−0.03	0.01	−0.01	9.48$\times {10}^{-3}$	9.15$\times {10}^{-3}$	9.46$\times {10}^{-3}$	9.42$\times {10}^{-3}$
	Figure 1viii	0.02	0.01	0.01	0.01	0.02	0.03	0.04	0.04	1.47$\times {10}^{-3}$	1.45$\times {10}^{-3}$	−4.02$\times {10}^{-3}$	4.10$\times {10}^{-3}$
$T_2^{†}A{T}_2$	Figure 1ix	1.19	1.19	1.48	1.47	1.66	1.66	1.84	1.85	0.82	0.82	0.98	0.97
	Figure 1x	−0.01	−0.01	0.00	0.00	0.00	0.00	0.00	0.00	9.92$\times {10}^{-4}$	1.57$\times {10}^{-3}$	−2.43$\times {10}^{-3}$	−2.99$\times {10}^{-3}$
	Figure 1xi	1.14	1.16	1.40	1.41	1.57	1.62	1.73	1.81	0.79	0.78	0.95	0.94
	Figure 1xii	−0.01	−0.01	0.00	0.00	0.00	0.00	0.00	0.00	9.22$\times {10}^{-4}$	1.57$\times {10}^{-3}$	−2.43$\times {10}^{-3}$	−2.99$\times {10}^{-3}$
	Figure 1xiii	−1.26	−1.25	−1.53	−1.52	−1.72	−1.70	−1.91	−1.89	−0.84	−0.84	−1.01	−1.01
	Figure 1xiv	0.01	0.01	0.01	0.01	0.00	0.00	0.01	0.00	7.35$\times {10}^{-3}$	7.37$\times {10}^{-3}$	−8.38$\times {10}^{-3}$	0.01
	Figure 1xv	−1.23	−1.21	−1.51	−1.48	−1.70	−1.67	−1.91	−1.85	−0.83	−0.83	−0.99	−0.98
	Figure 1xvi	0.01	0.01	0.01	0.01	0.00	0.00	0.01	0.00	7.35$\times {10}^{-3}$	7.37$\times {10}^{-3}$	−8.38$\times {10}^{-3}$	0.01
NC		771.15		929.91		1006.45		1079.44		378.74		581.00
Total		3.25	3.45	3.26	3.45	2.62	2.94	1.50	2.01	3.57	3.60	3.32	3.37

. The tables give the LERCCSD contributions, where the property vertex is O, as well as the nLERCCSD values, where the property vertex could be of the $p-p$, $h-h$, or the $p-h$ type, depending on the diagram. For example, the contribution from Sub-figure (ii) of Figure 1, for the nLERCCSD case, involved a p-h vertex, that is $O_{p-h}{T}_1$, and therefore included in it the contributions from the 26 diagrams in Figure 4. In general, O or $O_{x-y}$ (where “x” and “y” could be p or h) was the eEDM Hamiltonian for ${\mathcal{E}}_{\mathrm{eff}}$ (which is given in Table 2), while it was the dipole operator for the PDM (which is presented in Table 3).

Table 2 shows that for all the systems, the $A{T}_1$ term always dominated among the correlation terms, where A could correspond to either $O{T}_1$ or $O_{x-y}{T}_1$ for LERCCSD or nLERCCSD, respectively. For the effective electric fields of the HgX molecules, in the LERCCSD case, there were strong cancellations among the positive $A{T}_1$ and the negative $T_1^{†}A{T}_1$ and $T_2^{†}A{T}_2$ terms. However, the final values of nLERCCSD and the LERCCSD calculations matched within 2.5 percent, since in the former case, the $A{T}_1$ values were significantly lower than the latter, and the $T_1^{†}A{T}_1$ sector provided a far smaller contribution. In the case of SrF, the $A{T}_1$ terms were comparable for LERCCSD and the nLERCCSD scenarios, and therefore, inclusion of non-linear terms did not change its effective electric field, while for BaF, we observed a mechanism that was similar to that for the HgX systems. We observed a different trend for the PDMs, in Table 3. As the DF value and the nuclear contribution were the same for a given molecule, whether it be LERCCSD or an nLERCCSD calculation, the interplay between $A{T}_1$ and $T_1^{†}A{T}_1$ terms decided the importance of non-linear terms. For HgX, the $A{T}_1$ term in nLERCCSD calculations was always slightly larger in magnitude than the LERCCSD ones, while the net contributions from the $T_1^{†}A{T}_1$ terms, which were less significant, were the other way round. The resulting non-linear effects were not so important for SrF and BaF, as seen in the earlier paragraph, while for HgX molecules, it became significant, with their effects changing the PDM by up to about 25 percent for HgI.

We now conduct a survey of previous works on the effective electric fields and PDMs of the molecules that we considered, in Table 1. For the effective electric fields of BaF, we observed that the effective core potential-restricted active space SCF(ECP-RASSCF) [51] and restricted active space CI(RASCI) [52] methods gave larger values, while the result from the MRCIapproach in [54] estimated the values as being slightly lower, with respect to our nLERCCSD value. A discussion of the previous works on the effective electric fields of HgX and our improved estimate of the quantity using LERCCSD approach was already presented in [8], and hence, we re-direct the reader to the earlier work. Our nLERCCSD results improved over the earlier LERCCSD and FFRCCSD results, as both of those were calculated using a DZ quality basis. Most of the works that calculated PDMs and that were not mentioned in the table, including [39,41,59,60,61], were expounded in our earlier works in detail [9,33], and therefore, we only discuss in this paragraph the more recent works. The differences in the PDMs between the LERCCSD results in our earlier work and those in the present work for HgX were due to the choice of basis (DZ basis functions [8,57] in the former, as against a QZ basis in the current work). We observed that the values of PDM for SrF that were obtained by using the complete active space self-consistent field (CASSCF) approach to multi-reference CI (MRCI) and second-order Rayleigh–Schrodinger perturbation theory (RSPT2) [46] (which agreed with our nLERCCSD, as well as FFRCCSD results) underestimated and overestimated the results with respect to the experiment, respectively. The results for PDMs of SrF and BaF from Hao et al. [48] using the exact two-component Hamiltonian–Fock space coupled-cluster (X2C-FSCC) formalism and the PDM of SrF from Sasmal et al. [47] using a relativistic Z-vector coupled-cluster approach (with both works employing a QZ basis) agreed closely with the experimental values. However, we also note that while the Z-vector approach predicted the PDM of SrF very accurately, it underestimated that of BaF [50]. This existing difference in the PDMs of SrF and BaF between the nLERCCSD and the FFRCCSD approaches on one side and the Z-vector RCCSD approach on the other could possibly be resolved in future works that employ methods that are even more refined.

We now check for the scalability of our code that was parallelized using MPI. We did so by testing it with the SrF molecule, using a DZ basis. The code was to calculate both the effective electric field and PDM of the molecule for this test. As the code was structured in a way that the extent of parallelization was limited by the number of virtuals, which was 208 in this case, we chose to study scaling up to 192 processors (across eight nodes, and with 24 processors employed per node). The details of the computer (VIKRAM-100 super-computing facility at Physical Research Laboratory, Ahmedabad, India) that we used are: a 100-teraflop IBM nx360 M5 machine with 1848 processors. Each node had 24 processors (two Intel Xeon E5-2670 v3, each with 12 cores) and a memory of 256 GB RAM. The inter-process communication was via a 100% non-blocking FAT Tree Topology FDR (56 Gbits/s) Infiniband network. We used an Intel 15.2 compiler and impi5.0 and mkl libraries. As Figure 5 shows, our calculations indicated that the code was scalable up to this mark. In the figure, we plot the speedup against the number of processors, where the former was defined as $S_p=t_1/t_p$, with $t_p$ referring to the time taken for a computation with p processors. Performing the computations in serial mode took about 6.5 days, while calculations with four processors consumed around two days. The code took only 2.17 h to finish with 192 processors. The wall time approached saturation after 96 processors (2.51 h to 2.17 h from 96 to 192 processors), and hence, the optimal number of processors to use was around 96, where we still obtained a speed-up from 6.5 days to 2.51 h. The wall times were reliable as estimates, but not extremely accurate, as the computations were performed on a common cluster, and the speeds depended on other factors such as the number of users, the computer’s specifications, and type of jobs during the time interval across which our computations were done, although we took utmost care to ensure that no other application ran on the same node(s) as ours. However, our analysis was sufficient for the purposes of broadly demonstrating that our code was scalable to a reasonably large number of processors.

We also estimated the errors in our calculations. We first examined the error due to the choice of basis. We used QZ quality basis sets for our calculations, and as there was no five-zeta basis that was available for us to carry out any kind of estimate, we calculated the effective electric fields and the PDMs at the DZ level of the basis with our nLERCCSD code. We found that the percentage fraction difference between the DZ and QZ basis for ${\mathcal{E}}_{\mathrm{eff}}$ was around 3, 4, 5, and 7 percent for HgF, HgCl, HgBr, and HgI, respectively. We did not anticipate the difference between the DZ and QZ estimates to be over 10 percent for SrF and BaF either. Therefore, we did not expect that the difference between results from a higher quality basis set than QZ and those from a QZ basis should exceed 10 percent. Based on similar considerations, we estimated the error due to the choice of basis for the PDM to be at most 15 percent. Next, we shall look at the errors due to the ignored non-linear terms. They were expected to be negligible, and we ascribed a conservative estimate of two percent, which was the percentage fraction difference between the DF values and the current nLERCCSD values for the HgX molecules. Lastly, we comment on the importance of triple and other higher excitations. Based on our previous works and error estimates in them, we expected that these excitations would be around three percent for the purposes of calculating ${\mathcal{E}}_{\mathrm{eff}}$ [9], but could become important for PDMs. In conclusion, we linearly added the uncertainties and set an optimistic error estimate for the effective electric fields at about 15 percent. However, it was not so straightforward to set an error estimate for PDMs, as seen above, but we did not anticipate it to exceed 20 percent.

4. Conclusions

We investigated the contributions from the non-linear terms of the property evaluating the expression of the relativistic coupled-cluster theory in the determination of permanent electric dipole moments and effective electric field due to the electron electric dipole moment of SrF, BaF, and mercury monohalide (HgX with X = F, Cl, Br, and I) molecules. We found that the inclusion of these terms at the single and double excitation approximation brought the permanent electric dipole moments (PDMs) of SrF and BaF closer to the previously calculated finite-field relativistic coupled-cluster values, which were found to have overestimated the PDMs of the two molecules with respect to available measurements. The non-linear terms considerably changed the PDMs of HgX systems. For all of the chosen molecules, the non-linear terms were found not to change significantly the values of effective electric fields with respect to the results from the linear expectation value approach. However, such a result was a consequence of several cancellations at work. Since accurate estimation of these quantities are of immense interest to probe new physics from the electron electric dipole moment studies using molecules, our analysis demonstrated the importance of considering non-linear terms in relativistic coupled-cluster theory for their evaluations. We also presented the scaling behavior of our code with a representative SrF molecule and discussed the error estimates.