Relativistic Fock Space Coupled Cluster Method for Many-Electron Systems: Non-Perturbative Account for Connected Triple Excitations

The Fock space relativistic coupled cluster method (FS-RCC) is one of the most promising tools of electronic structure modeling for atomic and molecular systems containing heavy nuclei. Until recently, capabilities of the FS-RCC method were severely restricted by the fact that only single and double excitations in the exponential parametrization of the wave operator were considered. We report the design and the first computer implementation of FS-RCC schemes with full and simplified non-perturbative account for triple excitations in the cluster operator. Numerical stability of the new computational scheme and thus its applicability to a wide variety of molecular electronic states is ensured using the dynamic shift technique combined with the extrapolation to zero-shift limit. Pilot applications to atomic (Tl, Pb) and molecular (TlH) systems reported in the paper indicate that the breakthrough in accuracy and predictive power of the electronic structure calculations for heavy-element compounds can be achieved. Moreover, the described approach can provide a firm basis for high-precision modeling of heavy molecular systems with several open shells, including actinide compounds.


Introduction
The ultimate goal of modern relativistic quantum chemistry is the development of highly accurate ab initio methods of electronic structure modeling for studying atomic and molecular electronic structure, being applicable to systems with general geometry and shell structure, with emphasis on species containing heavy or super-heavy elements. In such systems relativistic and electron correlation effects are often not only of the same order of magnitude, but also non-additive and substantially intertwined. Systems which also possess strong quasi-degenerate or multi-reference character, i.e., cannot be described approximately by a single Slater determinant, are considerably more difficult for modeling, and conventional single-reference methods implemented in most common off-the-shelf packages are not always applicable. Some well-known atomic examples are Ni, Pd and Pt atoms, where the d 8 s 2 , d 9 s and d 10 electronic configurations are very close in energy to each other, and their intermixing presents serious problems to computational approaches and thus requires special treatment. A notable molecular example is U 2 , for which standard methods give unsatisfactory accuracy [1]. The first-principles-based methods aimed at high-precision modeling of such systems should not only ensure a balanced treatment of the dynamic and non-dynamic (e.g., caused by quasi-degeneracies) correlation effects, but also include relativistic and in many cases also Breit interaction and quantum electrodynamic (QED) effects on equal footing from the outset simultaneously to high order. Special relativity is fundamental for understanding chemical behavior of heavy elements [2,3], and the non-relativistic Hamiltonian can no longer serve as a good model for description of electronic subsystem and has to be replaced by other models derived from the relativistic consideration of particle motion, e.g., those based on the Dirac-Coulomb or Dirac-Coulomb-Breit Hamiltonians [4]. The most natural way to treat a manifold of quasi-degenerate problems such as excited states and bond breaking one should use multi-reference (MR) approaches, providing the flexibility needed to describe wavefunctions in such challenging cases [5]. For this reason, considerable efforts have been made to formulate MR methods and provide their implementations as general-purpose computer codes.
Ideally, these approaches should satisfy most or all of the following requirements: (i) size extensivity and size consistency; (ii) ability to handle various types of electronic shells; (iii) numerical stability in a wide range of molecular geometries; (iv) balanced treatment of all electronic states under consideration and predictable accuracy in all regions of potential energy surfaces; (v) acceptable computational cost. In addition the development of high-precision electronic structure methods is directed towards (i) improving the approximation of the Hamiltonian by inclusion of high-order relativistic and QED effects; (ii) improving the dynamic correlation treatment by increasing the excitation rank of cluster operators (in case of coupled cluster methods) or the perturbation theory (PT) order; and/or (iii) improving the non-dynamic correlation treatment by increasing size of the active valence or model (P) space, (iv) converging basis set to completeness.
Extensive studies over the past few decades have demonstrated that one of the most accurate (and thus appropriate for high-precision modeling of quasi-degenerate systems) method is the state-of-the-art "all orders" multi-reference CC (MRCC) method [6][7][8]. Some formulations of this method indeed meet all the requirements listed above. Genuine MRCC approaches can be divided into two main types, namely Fock space (FS) [9,10] and Hilbert space (HS) [11][12][13][14][15][16] approaches. Both types of MRCC formulations are based on the concept of multidimensional model space P, in which the effective Hamiltonian H eff is to be found by solving the generalized Bloch equation and then diagonalized [17,18]. Electron correlation is treated in two stages: dynamic correlation is described by virtual excitations to complementary Q-space determinants via the wave or cluster operator, while non-dynamic correlation is taken into account at the step of the H eff diagonalization. The distinction between these two kinds of correlation is rather arbitrary. The effective Hamiltonian MRCC approaches work very well in quasi-degenerate cases, but normally suffer from the convergence problems due to single or multiple intruder states. The intermediate Hamiltonian (IH) generalization of the effective Hamiltonian concept, which could be regarded as a good solution of the intruder state problem, was initially introduced by Malrieu et al. in the framework of degenerate PT [19]. This technique is very general, efficient and dramatically increases accuracy and area of applicability of the MR approach. Thus, the IH often becomes the method of choice for many challenging applications in the field of modeling of complicated atomic and molecular spectra within the both Fock space [20][21][22][23][24][25][26][27][28][29] and Hilbert space [15,16,30,31] MRCC approaches. Another efficient and rather universal way of avoiding the intruder state problem, related in a sense to the IH technique and perfectly fitting the FS-RCC framework, consists of using dynamic denominator shifts and extrapolated projected effective Hamiltonians [32,33].
Exploration of the multi-reference Hilbert space CC, which is generally more flexible than the FS-CC approach in the composing of model spaces for complex open shell electronic configurations, is still being rather limited due to the frequently occurring symmetry breaking problems (demanding additional non-trivial technical handling). The advantages of the Fock space MRCC method which led to its much more frequent use for solving real practical problems have roots in the fact that all excitation operators are defined with respect to the common vacuum vector. Unlike all HS MRCC approaches based on the Jeziorski-Monkhorst Ansatz, FS MRCC by construction do not suffer from the problem of redundancy of cluster amplitudes and hence preserves transparent physical picture intrinsic for the single-reference coupled cluster method as well as relative simplicity of computer implementation [34,35].
During recent decades the relativistic Fock Space CC approach in both effective and intermediate Hamiltonian formulations proved its efficiency and high accuracy in challenging calculations of heavy atomic and molecular systems with dense spectra, especially in cases when both d and f atomic orbitals become valence (see also recent reviews [26,29,36]). In the non-relativistic realm the pilot implementation of the FS-CC method accounting for triple excitations and allowing the treatment of up to six valence electrons sectors was pioneered by Hughes and Kaldor in earlier 1990s; they reported calculations of electronic states of light atoms in rather modest basis sets [37][38][39]. Later the more advanced non-relativistic implementation of the IH-based full FS-CCSDT approach (and some its effective approximations) to excitation spectra of some light atomic and molecular systems has been published by Musiał and co-authors [24,25,40]. However, to the authors' knowledge till now all the relativistic implementations of the FS-CC approach are limited to single and double excitations (i.e., CCSD approximation) and to the Fock space sectors with no more than two valence quasiparticles. It is to be mentioned that unlike its non-relativistic counterpart the FS-RCC approach has become one of the most widely used tools for relativistic electronic structure modeling. The list of fundamental research assisted by the FS-RCC predictions includes first experiments on laser resonance ionization spectroscopy of short-lived radioactive atoms and molecules [41][42][43] and searches for the electron electric dipole moment and related problems [44][45][46]. However, as was mentioned above, the accuracy and predictive power of FS-RCC was limited by the implementation of only the singles and doubles (CCSD) model up to the date [34]. High-precision atomic FS-RCCSD calculations have shown that even in the case of very extended basis sets and model spaces errors of order up to several hundreds of cm −1 remain for some complicated atomic cases and can be attributed only to the neglecting of higher rank excitations in the cluster operator [47]. Thus, the breakthrough in accuracy of high-precision electronic structure modeling can be achieved by systematic development of models including these contributions. Similar to the case of well-established non-relativistic single-reference coupled cluster theory, the development and applications of such the CC models to real problems are nearly impossible without efficient algorithms and computer codes due to the considerably high time complexity. Furthermore, any schemes of accounting for triples in some approximate manner (including calculations with truncated basis sets and spinor spaces) are also highly desirable, at least at the current state of hardware capabilities.
In this paper, we present the formulation of the relativistic Fock space coupled cluster method with iterative triples and its first program implementation allowing both atomic and molecular FS-RCCSDT calculations. Theoretical considerations are discussed in Section 2, and the pilot applications are reported in Section 3. Please note that not only energy calculations were conducted and reported (Sections 3.1 and 3.2), but also potential benefits for high-precision calculations of properties are revealed on the example of static dipole polarizability (Section 3.3).

Theory
We start with the conventional formulation of the FS-RCC approach described in details elsewhere [9,48]. This formulation implies the use of complete model spaces defined unambiguously by the choice of "active" (variable-occupancy) spinors and the universal exponential wave operator Ω normal-ordered with respect to the closed-shell Fermi vacuum state: The cluster operator T is a linear combination of excitation operatorsê K , where K stands for the list of normal-ordered creation/annihilation operators forming the excitation.
Electronic states energies and model-space parts of the corresponding wavefunctions are obtained through constructing and diagonalizing the effective Hamiltoniañ where H is the chosen approximation for the many-electron relativistic Hamiltonian, the subscript Cl marks the closed part of an operator and the overbar denotes its connected part. The coefficients t K (cluster amplitudes) should satisfy the equations which are usually written in the form where H 0 denotes the diagonal part of the Fock operator for the Fermi vacuum state. The entities D K (energy denominators) are the negatives of the differences of H 0 eigenvalues associated with the excitations K.
It is convenient to partition the cluster operator T according to the number of valence holes (n h ) and valence particles (n p ) to be created/destroyed (i.e., related to (n h h,n p p) sectors of the Fock space): According to the "subsystem embedding condition" [18], in order to describe the electronic states in the (N h h,N p p) sector of the Fock space, one needs to determine only T (n h h,n p p) with n h ≤ N h and n p ≤ N p . Therefore the system of coupled Equation (4) is split into subsystems, which can be solved consecutively. The T (n h h,n p p) operator at the CCSDT level of approximation includes connected single, double and triple excitations, Please note that the T (0h,0p) operator is actually the same that is used in the conventional single-reference CC theory. Diagrams representing cluster operators for the (0h,0p), (0h,1p) and (0h,2p) Fock space sectors are shown on Figure 1. The corresponding diagrams for other sectors are obtained in a straightforward manner by "turning down" valence lines. Please note that throughout the paper Hugenholtz-Brandow (antisymmetrized) diagrams are used (details of the formalism can be found elsewhere [49]). The well-known computational bottleneck of all CC models fully accounting for the triple excitations arises from the presence of the diagrams with contractions over four particle indices. These diagrams cause the N 8 scaling for compact systems where no advantages can be taken from localization techniques [50] (N is the number of one-electron basis functions). These diagrams for the (0p, 0h), (0p, 1h), and (0p, 2h) Fock space sectors are presented on Figure 2. Numerous approximate techniques reducing the overall computational cost of the CCSDT model have been proposed during recent decades, both within the single-reference [51][52][53][54][55] and multi-reference CC frameworks [37,40,[56][57][58][59][60][61][62]. The most widely used approximations are CCSD(T) and CCSDT-1, the former is often referred as the "gold standard" of quantum chemistry [63]. Here we generally follow the consideration proposed by A. Haque, S. Hughes and U. Kaldor [37,56,57] for the non-relativistic FS-CC to define various FS-RCC models accounting for triples. The general idea is to avoid the costly N 8 diagrams; a way for reaching this goal can be provided by the simple analysis based on many-body perturbation theory (PT) arguments. Since the T 2 operator arises in the first order of PT, the simplest estimate of the T 3 operator, is at least of the second order in perturbation V (see Figure 3). According to the Equation (3) for the effective Hamiltonian, triples will contribute in energy at least in the third PT order. If we evaluate triples amplitudes using the Equation (7) with singles and doubles amplitudes taken from converged FS-CCSD calculation, we arrive at the FS-CCSD+T(3) approximation which resembles the conventional CCSD(T) approach. However, despite its conceptual simplicity and low computational cost, the FS-CCSD+T(3) approach was shown to yield unpredictable results with worse accuracy than FS-CCSD [58,60,64]. The earlier conclusions by U. Kaldor and co-workers concerning the satisfactory performance of the FS-CCSD+T(3) scheme [56,57] were drawn from the direct comparison of the FS-CC excitation energy estimates, calculated with very small basis sets and thus necessarily unreliable, with their experimental counterparts. A straightforward improvement of this scheme consists of accounting for contributions from triple excitations in the r.h.s. of Equation (4) for single and double excitations amplitudes. Triples amplitudes should be re-calculated at each step of the iterative solution, resulting in increased computational cost, though the time complexity of the algorithms employed is actually the same (N 7 ) as in the FS-CCSD+T (3) scheme. This approximation introduced in [51,52] for the single-reference CC theory is referred to as the CCSDT-1 model (more precisely, the CCSDT-1b model). Again, it was extended to the non-relativistic multi-reference domain by U. Kaldor and co-workers and has shown rather promising results [37]; however, due to limitations caused by the high computational cost the features of the model were not well understood. Two other CCSDT-n models which are quite popular in the conventional CC theory, namely CCSDT-2 and CCSDT-3 [53], have never been generalized to the Fock space CC case. However, the generalization is straightforward: new diagrams corresponding to Ω ≈ e T 2 (FS-CCSDT-2) or Ω ≈ e T 1 +T 2 (FS-CCSDT-3) arise in the triples equations [49]. No PT-based restriction on diagrams is imposed, but the N 8 terms are avoided by simply prohibiting the diagrams with the triples amplitudes in the triples equations. The approximate treatment of the folded terms in the FS-CCSDT-n models requires special attention. The list of the folded diagrams arising in the full FS-CCSDT model is presented on the Figure 4. For consistency reasons only the terms appearing in the second PT order are to be retained within the FS-CCSDT-1 model. All such terms are represented by the (b) diagram arising in the (0h,2p) sector. We shall call the corresponding model FS-CCSDT-1 . Furthermore, taking into account the choice of the "direct" terms used in the T 3 amplitude equations within the FS-CCSDT-1 approximation (see Figure 3), it seems quite logical to retain only the leading (undressed) part of interaction in Figure 4b. In the latter case we arrive at a somewhat different model which will be referred to as FS-CCSDT-1 (to be consistent with the earlier work of U. Kaldor and co-authors).
To formulate the Fock space CCSDT-2,3 models, we should not use the PT arguments in order to be consistent with the single-reference counterparts, and the list of diagrams contributing to the triples equations should not include those ones containing any contractions with triples. As a result, only the   [32,33] offers the possibility to extend model spaces without encountering numerical instabilities caused by intruder states [65]. Thus, a partial incorporation of the contributions from the excitations normally considered to be triples is possible even within the CCSD model. However, significant deviations of resulting molecular property values from the experimental ones [47] or from full configuration interaction calculations [25,40] remain, indicating that some important excitations are missed even in the cases when large model (and/or intermediate) spaces are employed. Indeed, in the case of the (0h,1p) sector the most important triples amplitudes correspond to the determinants with two core holes which are doubly excited with respect to the model space (so-called spectator amplitudes). In the frame of the FS-RCCSD model these determinants are generated only by (disconnected) products of cluster operator terms inherited from the vacuum sector, so that the effect of differential correlation in the presence of the one electron added to the vacuum state is neglected. Similar problems are most pronounced for hole sectors, e.g., (1h,0p) [66]. For particle sectors the problem becomes even more severe in the case of (0h,2p) sector, thus we can expect the contribution of triples to be much more significant there.
All the electronic structure models discussed above were implemented in the EXP-T program package [35,67].

Pilot Applications
Here the pilot applications of the developed coupled cluster models are reported. Due to the high computational complexity of the schemes described above we employed a two-stage approach to achieve high accuracy sufficient for the direct comparison of calculated values with the experimental ones. At the first stage, FS-RCCSD calculations were performed using large flexible basis sets (denoted as LB) guaranteeing the smallness of the basis set incompleteness errors compared to the contribution from the triple excitations to the property under study; the calculations accounting for triples using such large basis sets are not feasible to the moment. At the second stage, the FS-RCCSD and the FS-RCCSDT-n/CCSDT calculations were performed within somewhat smaller basis (SB) in order to estimate the triples contribution to the desired quantity (excitation energy or other property). Obviously, model spaces employed in the LB and SB calculations were chosen consistently. Final estimates for the total energies were obtained as This naturally leads to similar estimates of transition energies and other properties considered in the present paper. The use of combined schemes like Equation (8) is a common and successful practice that has been used for different applications before [68][69][70][71][72]. High-order contributions calculated in such a way were shown earlier to be quite stable with respect to the choice of basis set and approximation to the relativistic Hamiltonian [69]. The functions of "small" basis sets were normally constructed (if not stated otherwise) from those of "large" ones as averaged atomic natural orbitals [73,74] obtained in scalar-relativistic CCSD calculations performed with the CFOUR program package [75]. Further details concerning the basis sets employed are provided in the corresponding subsections below; all basis set parameters are available online [76]. "Small-core" shape-consistent semilocal relativistic pseudopotentials by N. S. Mosyagin et al. [76][77][78] leaving both valence and sub-valence electrons of heavy atoms were used throughout. All explicitly treated electrons of heavy atoms as well as all electrons of light atoms were correlated.
Molecular integrals transformed to the basis of molecular spinors were calculated within the DIRAC package [79] and then imported to EXP-T.

Atomic Energy Levels of Thallium and Lead
High-precision calculations of atomic energy levels are of great importance for modern physics of super-heavy elements, e.g., for future spectroscopic experiments on atoms and elucidation of periodicity in their physical and possible chemical properties via revealing the nature of their ground and low-lying excited electronic states. To assess the accuracy of the model employed, calculations on lighter homologues are typically performed [47,80]. In particular, it was found that the errors in the FS-RCCSD predictions of atomic transition energies for thallium (∼80 cm −1 [80]) and lead (∼200-600 cm −1 [47]) can be attributed only to the neglect of higher excitations. Similar level of accuracy was also demonstrated recently by the MR CCSD-like approach of Dzuba and Flambaum [81]. We performed the series of FS-RCC calculations aimed at clarification of the role of higher excitations in the aforementioned problems.
High symmetry of atomic systems makes feasible calculations with very extended basis sets. Thus, the basis set limit in FS-RCCSD calculation can be nearly reached and the remaining error can be clearly attributed to the neglect of higher excitations. The FS-RCCSD calculations with the Dirac-Coulomb-Breit Hamiltonian were performed using the Tel Aviv Relativistic Atomic Fock Space coupled cluster code (TRAFS-3C) based on the radial symmetry, written by E. Eliav, U. Kaldor and Y. Ishikawa. The uncontracted universal basis set [82] used as the LB basis consists of even-tempered Gaussian type orbitals, with exponential parameters ζ given by the geometric series ζ n = αβ n−1 , α = 106, 111, 395.371615, β = 0.486752256286. The basis set used for both elements consists of 37s (n = 1 − 37), 32p (n = 5 − 36), 24d (n = 13 − 36), 21 f (n = 15 − 35), 12g (n = 21 − 32), 10h (n = 22 − 31), and 9i (n = 22 − 30) functions. All electrons were correlated, and virtual orbitals with energies above 300 a.u. were omitted. To account for the QED corrections to the energies we applied the model Lamb shift operator (MLSO) of Shabaev and co-workers [83]. The contributions from triples amplitudes to ionization potentials and excitation energies were calculated within the EXP-T program. Smaller contracted Gaussian basis sets (SB) [6s7p5d5 f 4g3h2i] (Tl) and [6s6p5d5 f 4g3h2i] (Pb) obtained by the ANO-type contraction procedure outlined above were employed for triples calculations; 60 core electrons of Tl and Pb were substituted by relativistic pseudopotentials as was described above.
The closed-shell 6s 2 configurations states of Tl + and Pb 2+ were taken as the reference for the Fock space sequences Tl + (0h, 0p) → Tl 0 (0h, 1p) (9) We restricted our attention to the electronic states described by the 6s 2 6p 1 (Tl) and 6s 2 6p 2 (Pb) leading configurations, thus the active spaces in both LB and SB calculations comprised 6p-orbitals only.
The ionization potentials (IP) and low excitation energies (EE) of Tl, Pb and the Pb + cation are reported and compared with the experimental data in Table 1. Table 1. Deviations of the calculated ionization potentials (IP) and excitation energies (EE) of neutral thallium and lead and lead cation (cm −1 ) from the experimental values. FS-RCCSD/LB+T/SB stands for the combined scheme (8). One can state that even combined schemes with FS-RCCSDT-1 and FS-RCCSDT with minimal model spaces greatly outperforms the IH-FS-RCCSD model with large main and intermediate model spaces used for Pb [47]. This fact confirms the impossibility to recover the effects of all triples via extending model spaces in FS-RCCSD calculations. Moreover, in effectively single-reference problems like thallium atom model space enlargement does not lead to any improvement.

State Exptl IH-FS-FS-FS-RCCSD/LB + T/SB [84] RCCSD [47] RCCSD/LB SDT-1 SDT-1' SDT-2 SDT-3 SDT
The largest deviation from the experimental excitation energy for the FS-RCCSDT model was obtained for the rather high-lying 1 S 0 state of the Pb atom (174 cm −1 ). This can be explained by a strong interaction with the 6s 2 6p7p configurations which were not included into the model space in current calculations. The remaining systematic errors for other electronic states are less than 30 cm −1 and can be caused by imperfection of the additive scheme (8) used and by the neglect of quadruple and higher excitations in the cluster operator. However, the present results are to our knowledge the most accurate ab initio data for heavy non-alkali atoms. Thus, the technique described seems to provide the basis for the long-awaited breakthrough in accuracy at least in the domain of atomic electronic structure calculations.

Electronic States of TlH
The TlH molecule provides a good test case for the performance of relativistic electronic structure methods, because even qualitative description of its electronic states is impossible without careful treatment of the spin-orbit coupling. For this reason, the TlH molecule was often used for the benchmarking of new electronic structure models (see [85] and references therein). Low-lying electronic states of TlH converge to the Tl( 2 P o 1/2 ) + H( 2 S 1/2 ) and Tl( 2 P o 3/2 ) + H( 2 S 1/2 ) dissociation limits, and two of them (B1 and (1)0 − ) corresponding to the former limit are unstable. Reliable experimental data are available for the A0 + state [86,87].
Due to the relativistic stabilization of the 6s shell of thallium it seems natural to describe the electronic states of TlH in a wide range of internuclear separations in the (0h,2p) FS-RCC sector, considering the TlH 2+ ion as the Fermi vacuum. However, the optimal choice of active set of spinors is quite complicated. The most natural choice at large internuclear distances is the active space comprised of 4 lowest lying Kramers pairs of virtual spinors, corresponding to the 6p and 1s spinors of Tl and H atoms, respectively. However, at the internuclear distance comparable with the bond lengths (R ∼ 1.9 Å) the fourth pair of spinors becomes purely antibonding and appears very close in energy to higher one-electron states. As a result, no stable physically reasonable solutions of the FS-RCC equations can be obtained for this range of R. Another possible way is to include additional virtual spinors to the active space. However, the one-electron spectrum of TlH 2+ is very dense, with a complex dependence on the internuclear separation (numerous crossings and avoided crossings). To ensure a regular behavior of computed potential curves, one would need to include several dozens of Kramers pairs into the active space and, hence, relatively large values of denominator shifts [32] would be required to achieve stable solutions. This strategy leads to quite a blurry physical picture and to an unavoidable increase of extrapolation error for large shift amplitudes. Specifically for this problem the simplest solution is to restrict attention only to the "bonding" region and to choose only the 3 lowest virtual Kramers pairs as active; contributions of the determinants missed in the case of such a modest active space can be retained within the FS-RCCSDT model.
The scheme (8) was again employed to account for connected triple excitations. The [6s7p5d5 f 4g3h2i] (Tl) and aug-cc-pVTZ [88] (H) contracted basis sets were used in the SB calculations, while the uncontracted basis set (14s13p10d7 f 6g5h5i) (Tl) and the uncontracted version (7s4p3d2 f ) of the aug-cc-pVQZ set were used in the large basis (LB) FS-RCCSD calculations. 60 core electrons of Tl were represented by the relativistic pseudopotential (see above). During the calculations accounting for triples excitations the 5s and 5p shells of Tl were excluded from the correlation treatment.
To suppress the instabilities in the FS-RCC equations the small amplitudes of denominator shifts [32] in the non-trivial FS sectors were used, 0.1, 0.2 and 0.3 a.u. for singles, doubles and triples amplitudes, respectively; effective Hamiltonians computed at the values of the attenuation parameter n = 2-4 (see Equations (3) and (4) in Reference [33]) were extrapolated to the zero-shift limit by the [0/1] matrix Padé approximant [33]. The basis set superposition error was estimated by the counterpoise correction. Equilibrium distances as well as harmonic vibrational frequencies were evaluated with the help of the program VIBROT [89]. Equilibrium distance r e for the A0 + state was extracted from the rotational constants reported in [87] for TlD. Table 2 summarizes term energies and spectroscopic constants for the X0 + and A0 + states of the TlH molecule. Term energy of the A0 + state as well as equilibrium internuclear distances obtained by the complete ab initio and the triples-corrected models are in reasonable agreement with the experimental ones, while those derived from the FS-RCCSD potential curves can hardly be considered to be satisfactory. One can see that the inclusion of connected triples amplitudes brings the vibrational constants to a much better agreement with the experiment for both electronic states.
Further refinement of the term energy and spectroscopic constants of the A0 + state is possible by combining the ab initio transition energies as functions of the internuclear separation with accurate Rydberg-Klein-Rees (RKR) ground-state potential energy function derived from the spectroscopic constants reported in [90] (the RKR1 program by LeRoy [91] was employed). Resulting values obtained in the FS-RCCSDT-corrected calculation are in excellent agreement with the experimental ones; note that this result was obtained for the smallest possible model space for which the FS-RCCSD model yields very poor results.
Further experience with more complicated molecules is required to turn FS-RCCSDT into a standard tool in high-precision molecular calculations.

Static Dipole Polarizability of Lead
The developed method allows the calculation of different properties of atoms and molecules using the finite-field approach. One of such properties is the static dipole polarizability. It characterizes the interaction of an atom with an external electric field. It can be also used to estimate some chemical properties of an element such as its dispersion interactions with surfaces.
The polarizability tensor is defined as a mixed derivative of the energy with respect to components of the external electric field F via the following expression: We focus here on the isotropic (scalar) part of the tensor. Polarizability can be calculated numerically by the direct use of the finite-field procedure based on the Equation (11). For this purpose, one can turn on the interaction of the atom with the external electric field adding this interaction as a perturbation to the electronic Hamiltonian. In the case of an atom, such perturbation reduces the symmetry of the Hamiltonian from the spherical to the axial one. Since it is not necessary to employ the spherical symmetry in the calculations carried out with the EXP-T code, it is possible to avoid using any special procedure, such as the sum-over-states approach. Other properties such as the electron electric dipole moment enhancement factor can be calculated in a similar way [72,92,93]. Table 3 presents the calculated values of the polarizabilities of the ground electronic state of the lead atom and its ions. The leading contribution to the polarizability has been calculated at the FS-RCCSD level with the uncontracted basis set (LB) consisting of (25s25p22d10 f 6g5h4i) functions with exponential parameters forming the geometric progression. The use of such a large basis set is important to ensure the satisfactory treatment of all leading polarization effects. To calculate the contribution of triple cluster amplitudes (see Equation (8)) the SB basis set consisting of (13s11p10d2 f 2g) functions has been employed; the f and g functions were again contracted using the ANO-type procedure (see above). In all polarizability calculations 60 core electrons of Pb were substituted by the relativistic pseudopotential [77]. Table 3. Static scalar dipole polarizabilities for lead and its ions (in atomic units).

FS-FS-RCCSD/LB + T/SB RCCSD/LB SDT-1 SDT-1' SDT-2 SDT-3 SDT Experiment
Pb 2+ 6s 2  One can see from Table 3 that triples are mostly important for the neutral lead atom. In this case, triples contribution is about 7%. It can be noted that the values of the polarizability obtained within the FS-RCCSDT-3 model are very close to those obtained within the full FS-RCCSDT method.
The calculated values of the polarizability for different charge states of Pb are in good agreement with available experimental data (see Table 3). The polarizability of the ground state of the neutral lead is also in a reasonable agreement with previous theoretical data: 46.5 a.u. [98], 44.04 a.u. [81], 47.70 [96], 46.96 [99] (see Ref. [97] for a complete overview of the polarizability calculations). The polarizability of the Pb + cation has not ever been measured experimentally, but was theoretically estimated earlier to be equal to 23.5 a.u. [100]. We do not know about any other calculation of polarizabilities of Pb + and Pb 2+ .

Concluding Remarks
Pilot applications of the presented FS-RCCSDT and FS-RCCSDT-n schemes to atomic (Tl, Pb) and molecular (TlH) systems clearly demonstrate the importance of accounting for triple excitations in relativistic coupled cluster models for achieving really high accuracy. The results indicate that the decisive breakthrough in the accuracy and predictive power of the electronic structure modeling of heavy-element compounds can be achieved by the proposed way. However, further developments in both the theoretical foundations and the algorithm design are clearly required to make the FS-RCCSDT models applicable to wide classes of complex heavy-atom molecules as well as to cluster models of solids [101].
Finally, we note that the proposed models with non-perturbative triples provide a firm basis for the proper extension of the FS-RCC method to higher Fock space sectors (n h + n p ≥ 3), for instance, to the (1h,2p) and (0h,3p) sectors [37,39,102]. The restriction of the cluster expansions to singles and doubles in these cases actually means the building of the effective Hamiltonian for the target sector exclusively from the lower sectors diagrams with no more than two quasiparticles, and thus resulting in a low precision in most cases. This conclusion is supported by the numerical evidence reported in Reference [103], which clearly demonstrated the insufficiency of the triple-electron-affinity equation-of-motion (TEA-EOM) CCSD model for a quantitative description of systems with three open shells. The extension of the FS-RCCSDT model to higher FS sectors is under testing now and will be reported in our future papers.