Significance of non-linear terms in the relativistic coupled-cluster theory in the determination of molecular properties

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 singles and doubles excitations (RCCSD method). The effective electric fields (${\cal 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 ${\cal E}_{eff}$ and PDMs for the above molecules, we had 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 above quantities and their computational time scalability with 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.


I. 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 Refs. [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. Upper bound 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]. A 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, E 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, Ref. [5]), or to propose a new candidate for molecular eEDM experiments (for example, Ref. [8]). This quantity can only be obtained using a relativistic many-body theory [29]. Calculating the PDM provides information on polarizing factor for molecules that are proposed for an eEDM searches, where the property has not been measured.
as LERCCSD method). 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 1 percent) in the values of E eff . The results for the PDMs obtained in these methods were comparable for SrF and BaF and also were over-estimating 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 eventually be compensated for with further iterations.
Here, we intend to calculate the values of E ef f 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 Refs. [36,37] to implement these non-linear RCC terms. For this purpose, we have undertaken molecules that are 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 proposes to laser-cool HgF has opened new avenues for an upcoming eEDM experiment with the molecule [38]. Another very important molecule in this regard is BaF, and two eEDM experiments are simultaneously underway for this system [39,40]. Experimental values of the PDMs are available only for BaF among the systems that we mention 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 [18].

II. THEORY AND IMPLEMENTATION
In the RCC theory, the wave function of a molecular state is expressed as [41] |Ψ = e T |Φ 0 , where T is the cluster operator and |Φ 0 is the reference state obtained by mean-field theory. We use the Dirac-Coulomb Hamiltonian in our calculations, and |Φ 0 is obtained using the DF method. In the RCCSD method, we approximate T = T 1 + T 2 with subscripts 1 and 2 indicating singles and doubles excitations, respectively, and they FIG. II: The effective one-body terms representing particle-particle (p-p) diagrams considered in this work. i, j, k, · · · and a, b, c, · · · refer to holes and particles, respectively. The symbol of the operator, O p−p , is not mentioned explicitly in the diagrams, and the property vertex is the dashed line ending with an 'o' in each diagram. are given using the second-quantization operators as where the notation i, j is used to denote holes a, b refer to particles, t a i is the one hole-one particle excitation amplitude and t ab ij is the two-hole two-particle excitation amplitude. We have employed the UTChem [42,43] program for the DF calculations, the atomic orbital to molecular orbital integral transformations as well as generating the property integrals, and the Dirac08 [44] 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 where the subscript, 'c', means that each term is fully contracted [45], or in the diagrammatic terminology, connected [37]. The PDM of a molecule is determined as [16] where D is the electric dipole operator, the index A runs over the number of nuclei, Z A is the atomic number of the A th nucleus, and r A is the position vector from the origin to the site of the A th nucleus. The first term in the above expression is the electronic term, while the second term is the nuclear contribution. Similarly, E eff is evaluated as where the summation is over the number of electrons, N e , β is one of the Dirac matrices (also known as γ 0 in literature), Σ is the (4x4) version of Pauli matrices, and E intl i is the internal electric field that is experienced by the i th 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, 1 rij , and is complicated, we take recourse to employing an effective eEDM Hamiltonian instead of the one introduced above [29]. It follows that where γ 5 is the product of the gamma matrices (given by iγ 0 γ 1 γ 2 γ 3 ), while p i is the momentum of the i th electron.
In the LERCCSD method, the following expression has been used in the evaluation of the expectation values These terms are represented using Goldstone diagrams and are shown in Fig. I. Note that OT 2 and its hermitian conjugate are zero, due to Slater-Condon rules [46,47]. Diagrammatically, such a diagram will have at least two open lines, that is, it is not fully connected. The evaluation of properties using the LERCCSD approximation misses 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 Eq. (3) even in the RCCSD method as it contains a non-terminating expression. However, it is possible to demonstrate the importance FIG. IV: The list of the effective one-body terms representing the particle-hole (p-h) diagrams in this work. The notations are the same as in the particle-particle and the hole-hole diagrams.
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 Eq. (3) in heavier molecules, due to amount of computations involved in it. In order to tackle this issue, we adopt an additional computational step by breaking the non-linear terms into intermediate parts as described more elaborately below. Further, we parallelize 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 Fig. I. 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) too. 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 Fig. IV. We assign the notation O p−h for such a vertex. This sequence of operations already generates 26 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 hermitian conjugate of sub-figure (i) of Fig IV gives  As can be seen from the above discussions, some of the diagrams that are undertaken in this procedure demand up to the order of . This becomes especially relevant when we perform computations on heavy systems and with high-quality basis sets, such as those that we have chosen 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 is a full 5 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 can be seen from the heaviest HgI too. We add at this point that we have exploited the C 8 double point group symmetry in our nLERCCSD code, as we had done for the earlier LERCCSD program [5,48]. This aspect also substantially lessens the computational cost, as it restricts the number of matrix elements to be evaluated based on group theoretic considerations. For example, OT 1 involves computing matrix elements of the form a|O|i . Given that we have 89 holes and 429 particles, the number of possible matrix elements are ∼ 3.8 × 10 5 . However, since we impose the restriction that both i and a should belong to the same irreducible representation, we need to evaluate only ∼ 7.2 × 10 4 matrix elements. Similar considerations for the more complicated terms involving T 2 leads to evaluating much fewer matrix elements.  [49][50][51][52][53][54]. It is to be noted that the chosen values for the HgX molecules are from theory, while those for SrF and BaF are from experiment. Also, we opted for Dyall's quadruple zeta (QZ) basis for Hg and I [55], Dunning's correlation consistent polarized valence quadruple zeta (cc-pVQZ) basis for the halide atoms (F, Cl, and Br) [56], and Dyall's QZ functions augmented with Sapporo's diffuse functions [57] for Sr and Ba. We chose Dyall's basis for Hg and I as it is among the most reliable and widely used basis functions for heavy atoms. We did not add diffuse functions as it increases the computational cost drastically for QZ quality basis sets. Moreover, it was found that inclusion of diffuse functions change the effective electric field by around 2.5 percent for HgF, and is expected to lead to a similar difference for the heavier HgX systems [9]. However, in the foreseeable future, such computations could be peformed 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 can expect that the missing contributions would be minimal, and possibly even negligible.

III. RESULTS AND DISCUSSIONS
In Table I, we present our results for HgX, SrF, and BaF, all using QZ basis sets. We discuss the trends in the PDMs and E eff s across HgX in the three different approaches, namely LERCCSD, and nLERCCSD methods, while briefly making a comparison with the FFRCCSD method from Ref. [9] wherever relevant, and also examine the correlation effects in the property from lower to all-order methods. SrF and BaF molecules are treated as standalone systems. Firstly, we observe that the effect of non-linear corrections is to increase the PDM and decrease the effective electric field (except in the case of SrF, where the difference is still within 0.5 percent). We find from Table I that for SrF and BaF, the nLERCCSD method yields PDMs that are very close to their LERCCSD counterparts (within 1.5 percent of each other for both the molecules), but are not in better agreement with experiments than their LERCCSD counterparts. However, the nLERCCSD values agree 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 cannot be made with the HgX molecules, as available FFRCCSD data uses a double zeta (DZ) quality basis. For HgX systems, we observe that unlike in the cases of SrF and BaF, the difference between the LERCCSD and the nLERCCSD results widen from about 6 percent for HgF and HgCl, to about 25 percent for HgI. The values for E eff for SrF and BaF show that the LERCCSD, nLERCCSD, and FFRCCSD methods all agree to within 1 percent. In the case of HgX molecules, the LERCCSD and the nLERCCSD results are 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 found that its effective electric field is 110.87 GV/cm, which is lesser than the nLERCCSD value by 2.5 percent.
The individual contributions that arise from diagrams given in Fig. I to the effective electric fields and PDMs of HgX, SrF, and BaF molecules are given in Tables II and III Table II), while it is the dipole operator for the PDM (which is presented in Table III). Table II shows that for all the systems, the AT 1 term always dominates among the correlation terms, where A could correspond to either OT 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 are strong cancellations among the positive AT 1 and the negative FIG. V: Plot showing the scaling behaviour of the program in the property evaluating expression for a representative system, SrF, with number of processors of our computer. The X-axis gives the number of processors, while the Y-axis is the speedup, S p = t 1 /t p , where t is the time taken and the subscript denotes the number of processors. We have used a double-zeta quality basis for this purpose, and test up to 192 processors, as the parallelism in our code is limited by the number of virtual orbitals, which is 208 in this case.
T † 1 AT 1 and T † 2 AT 2 terms. However, the final values of nLERCCSD and the LERCCSD calculations match within 2.5 percent, since in the former case, the AT 1 values are significantly lower than the latter, and the T † 1 AT 1 sector provides a far smaller contribution. In the case of SrF, the AT 1 terms are comparable for LERCCSD and the nLERCCSD scenarios, and therefore inclusion of non-linear terms does not change its effective electric field, while for BaF, we observe a mechanism that is similar to that for the HgX systems. We observe a different trend for the PDMs, in Table III. As the DF value and the nuclear contribution are the same for a given molecule, whether it is LERCCSD or an nLERCCSD calculation, the interplay between AT 1 and T † 1 AT 1 terms decide the importance of non-linear terms. For HgX, the AT 1 term in nLERCCSD calculations is always slightly larger in magnitude than the LERCCSD ones, while the net contributions from the T † 1 AT 1 terms, which are less significant, are the other way round. The resulting non-linear effects are not so important for SrF and BaF as seen in the earlier paragraph, while for HgX molecules, it becomes 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 have considered, in Table I. For the effective electric fields of BaF, we observe that effective core potential-restricted active space SCF (ECP-RASSCF) [21] and restricted active space CI (RASCI) [22] methods give larger values, while the result from MRCI approach in Ref. [24] estimate 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 has already been presented in Ref. [8], and hence we re-direct the reader to the earlier work. Our nLERCCSD results improve over the earlier LERCCSD and FFRCCSD results, as both of those were calculated using a DZ quality basis. Most of the works that calculate PDMs and which are not mentioned in the table, including Refs. [58][59][60][61][62], have been expounded in our earlier works in detail [9,16], 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 are due to the choice of basis (DZ basis functions [8,27] in the former, as against a QZ basis in the current work). We observe that the values of PDM for SrF that are obtained by using complete active space self consistent field (CASSCF) approach to multi-reference CI (MRCI) and second-order Rayleigh-Schrodinger perturbation theory (RSPT2) [14] (which agrees with our nLERCCSD as well as FFRCCSD results) underestimate and overestimate the results with respect to experiment, respectively. The results for PDMs of SrF and BaF from Hao et al [17] using exact two-component Hamiltonian-Fock Space coupled-cluster (X2C-FSCC) formalism and the PDM of SrF from Sasmal et al [15] using a relativistic Z-vector coupled-cluster approach (with both the works employing a QZ basis) agree closely with experimental values. However, we also note that while the Z-vector approach predicts the PDM of SrF very accurately, it underestimates that of BaF [20]. 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 do so by testing it with the SrF molecule, using a DZ basis. The code is to calculate both effective electric field and PDM of the molecule for this test. As the code is structured in a way that the extent of parallelization is limited by the number of virtuals, which is 208 in this case, we chose to study scaling up to 192 processors (across 8 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 has 24 processors (two Intel Xeon E5-2670 v3, each with 12 cores) and a memory of 256 GB RAM. The inter-process communication is via a 100% non-blocking FAT Tree Topology FDR (56 Gbits/Sec) Infiniband network. We use an Intel 15.2 compiler, and impi5.0 and mkl libraries. As Figure V shows, our calculations indicate that the code is scalable up to this mark. In the figure, we plot the speedup against the number of processors, where the former is 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 takes about 6.5 days, while calculations with 4 processors consumes around 2 days. The code takes only 2.17 hours to finish with 192 processors. The walltime approaches saturation after 96 processors (2.51 hours to 2.17 hours from 96 to 192 processors), and hence the optimal number of processors to use is around 96, where we still get speed-up from 6.5 days to 2.51 hours. The walltimes are reliable as estimates, but not extremely accurate, as the computations were performed on a common cluster, and the speeds depend upon other factors such as the number of users, the computer's specifications, and type of jobs during the time interval across which our computations are done, although we took utmost care to ensure that no other application ran on the same node(s) as ours. However, our analysis is sufficient for the purposes of broadly demonstrating that our code is scalable to a reasonably large number of processors.
We also estimate the errors in our calculations. We first examine the error due to choice of basis. We use QZ quality basis sets for our calculations, and as there is no 5-zeta basis that is available for us to carry out any kind of estimate, we calculate the effective electric fields and the PDMs at DZ level of basis with our nLERCCSD code. We find that the percentage fraction difference between DZ and QZ basis for E eff is around 3, 4, 5, and 7 percent for HgF, HgCl, HgBr, and HgI, respectively. We do not anticipate the difference between the DZ and QZ estimates to be over 10 percent for SrF and BaF either. Therefore, we do 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 estimate the error due to 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 are expected to be negligible, and we shall ascribe a conservative estimate of 2 percent, which is 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 expect that these excitations would be around 3 percent for the purposes of calculating E eff [9], but could become important for PDMs. In conclusion, we linearly add the uncertainties and set an optimistic error estimate for the effective electric fields at about 15 percent. However, it is not so straightforward to set an error estimate for PDMs, as seen above, but we do not anticipate it to exceed 20 percent.

IV. CONCLUSION
We have investigated contributions from the non-linear terms of the property evaluating expression of the relativistic coupled-cluster theory in the determination of permanent electric dipole moments and effective electric field due to electron electric dipole moment of SrF, BaF and mercury monohalides (HgX with X = F, Cl, Br, and I) molecules. We find that the inclusion of these terms at the singles and doubles excitations approximation brings the permanent electric dipole moments (PDMs) of SrF and BaF closer to the previously calculated finite-field relativistic coupledcluster values, which were found to have overestimated the PDMs of the two molecules with respect to available measurements. The non-linear terms considerably change the PDMs of HgX systems. For all of the chosen molecules, the non-linear terms are found to not significantly change the values of effective electric fields with respect to the results from the linear expectation value approach. However, such a result is 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 demonstrates importance of considering non-linear terms in relativistic coupled-cluster theory for their evaluations. We have also presented the scaling behaviour of our code with a representative SrF molecule, and discussed the error estimates.