Triple excitations in Green's function coupled cluster solver for studies of strongly correlated systems in the framework of self-energy embedding theory

Embedding theories became important approaches used for accurate calculations of both molecules and solids. In these theories, a small chosen subset of orbitals is treated with an accurate method, called an impurity solver, capable of describing higher correlation effects. Ideally, such a chosen fragment should contain multiple orbitals responsible for the chemical and physical behavior of the compound. Handing a large number of chosen orbitals presents a very significant challenge for the current generation of solvers used in the physics and chemistry community. Here, we develop a Green's function coupled cluster singles doubles and triples (GFCCSDT) solver that can be used for a quantitative description in both molecules and solids. This solver allows us to treat orbital spaces that are inaccessible to other accurate solvers. At the same time, GFCCSDT maintains high accuracy of the resulting self-energy. Moreover, in conjunction with the GFCCSD solver, it allows us to test the systematic convergence of computational studies. Developing the CC family of solvers paves the road to fully systematic Green's function embedding calculations in solids. In this paper, we focus on the investigation of GFCCSDT self-energies for a strongly correlated problem of SrMnO$_3$ solid. Subsequently, we apply this solver to solid MnO showing that an approximate variant of GFCCSDT is capable of yielding a high accuracy orbital resolved spectral function.

Realistic molecules and solids often display features allowing them to be described by embedding methods. In these realistic problems, the orbital space usually can be divided into two parts: a fragment (system) requiring an accurate description and an environment that can be represented only approximately. Usually, the embedding construction is done in such a way that first the entire problem is evaluated via a low-scaling method yielding a low-energy impurity Hamiltonian for the system/chosen fragment. Afterwards, the low-energy system Hamiltonian is treated with a more accurate non-perturbative method.
Such an embedding construction is not only applicable for solids and molecules containing transition metals with d-and f -shells, where the local impurity is small and composed of d-and f -orbital, but also for dispersion dominated solids, for example, molecular crystals, two-dimensional layered materials. In these materials, the local impurity may be large and may involve a large number of orbitals, often a big supercell [1].
The choice of the non-perturbative method employed in the treatment of the chosen fragment (or the impurity problem) is often challenging for both SEET and DMFT. This is because this method (frequently called a solver) * ashee@berkeley.edu should treat the impurity problem accurately. In particular, for strongly correlated problems, it is hardly possible to treat impurities that are larger than 16 orbitals.
One of the most widely employed solvers, continuoustime quantum Monte Carlo (CTQMC) [14], is applicable only for a few sites, limited to Hubbard-type interactions, and fraught with fermionic sign problem. Another method, Exact Digaonalization (ED) [15], which requires discretized impurity bath can treat generalized interaction, but also limited to very few sites because of its exponential scaling. In order to alleviate those limitations many truncated wave-function based solvers, for example, DMRG-based [16], selected and truncated CI-based [17,18] solvers have been proposed. Usually these truncated solvers are employed for those problems, where a larger number of sites is considered and the accuracy can be somewhat sacrificed in comparison to ED.
We have chosen the coupled cluster [19,20] (CC) hierarchy of methods as solvers for solving the impurity problems. CC methods are very successful in quantum chemistry [21]. Recently, we have implemented a Green's function mapping of the CC wave function for use in embedding frameworks [22]. In our previous works, we have employed coupled cluster singles doubles truncation (CCSD) as an impurity solver called here GFCCSD. Similar solver has been developed by Zhu et. al. [23,24] in the dynamical mean field theory (DMFT) context. However, for a truncated solver it is not guaranteed to be successful at all interaction strengths. Therefore, in our previous works, both with periodic solids [25] and molecules [26] we have explored strengths and weaknesses of the GFCCSD solver. We observed that GFCCSD (i) frequently failed to identify a correct ground state of the impurity problem (ii) was reasonably successful when the number of impurity orbitals was small (iii) with a limited truncation at the SD level is hard to be completely decisive if the convergence with respect to excitation level was reached.
Therefore, in this work we will improve upon the GFCCSD solver by including triple excitations. For CC hierarchy, the inclusion of higher rank excitations provides a much faster convergence in comparison to other wave function methods, such as, configuration interaction (CI) and many body perturbation theory (MBPT) [27]. Moreover, historically the inclusion of triple excitations in schemes such as CCSD(T) provided unprecedented accuracy in molecular problems. With the inclusion of full triple excitations, the computational complexity increases from O(N 6 ) to O(N 8 ). Therefore, it is important to make approximations so that the steep scaling of a triples calculation can be reduced. We have proposed such an approximate scheme in the next section that has the computational complexity O(N 7 ). To our knowledge, this is the first reported GFCCSDT implementation including full triples, which has been used in a quantum embedding context. A GFCC scheme including approximate triples excitation has been reported previously by Peng and Kowalski [28].
The GFCC solver with triple excitations can be essential for quantitatively accurate embedding calculations. This is because in realistic applications of embedding methods, we often want to achieve a systematic improvement/convergence of the observables (total energy, spectral function etc.) in terms of various computational choices. Two such choices relevant only for the embedding part of a calculation are: (a) the accuracy of the bath representation in the impurity problem ; (b) the number of orbitals in the local problem chosen for embedding.
The GFCCSD solver with triple excitations will allow us to employ a very large number of bath orbitals, thus ensuring accurate bath representation and a convergence with respect to the number of bath orbitals. Moreover, it is often difficult to choose a local fragment which is optimal for the description of the chosen system. In GFCCSDT, we can systematically increase the number of orbitals that are chosen while still maintaining high accuracy of calculations. Both of these advantages of the newly developed solver will allow for quantum embedding calculations to become systematically improvable and highly quantitatively accurate.
CC is a many-body theory based on the exponential parametrization of the ket wave function |Ψ = eT |Φ , where |Φ is a reference mean-field wavefunction. In this paper, we are using a unrestricted Hartree Fock (UHF) determinant |Φ UHF as the reference wave function. T is a cluster operator inducing various h-p excitation from the reference wave function |Φ .
In CC, the bra wave function is not an adjoint of the ket because eT is not a unitary operator. In general, the choice of the bra state is non-unique in CC. Here, we define Ψ L | = Φ|(1 +Λ)e −T through the de-excitation operator Λ as the bra state [29]. This definition ensures that the bra and ket states are biorthonormal. A CC wave function is mapped to a corresponding Green's function [22,[30][31][32] by the following definition: where a p = e −T a p eT and a † p = e −T a † p eT . The creation (annihilation) operator is a † p (a p ) operating on a singleparticle state p. The transformed Hamiltonian is defined as H = e −T HeT − E gr . The UCC ground state energy is denoted as E gr . While CC is equivalent to ED if no truncation is made in the rank of the cluster opera-torT , in practical implementations, the rank of cluster operator is always truncated to minimize computational cost. In this work, we employ the singles-doubles-triples (SDT) approximation, resulting in the following opera-torsT =T 1 +T 2 +T 3 ;Λ =Λ 1 +Λ 2 +Λ 3 for the ground state problem. In a subsequent step, we tridiagonalize H in the space of (N + 1) and (N − 1) electronic wave functions, using the Lanczos algorithm. The (N + 1) and (N − 1) electronic states are spanned by all the functions generated byR andL operators defined aŝ and acting on eT |Φ and Φ|e −T , respectively. BothL andR operators are truncated at the level of triples.
The computational scaling of the evaluation of the ground state problem with SDT approximation scales as n 8 b , whereas the tridiagonalization step scales as where n b is the number of basis functions. The maximum memory requirement for the SDT truncation is due to storing T 3 and Λ 3 amplitudes. In order to make both computational scaling and memory requirements manageable, methods that include triples are often employed with various approximations that effectively lower the compuational and memory scaling. In this work, we have used an approximation proposed by Hirata et. al. [33] in their EOM-CCSD (m,n) work, where m and n represents the rank of excitations in the ground state and in the excited state problem, respectively. Krylov and Slipchenko [34] later analyzed the strength and weaknesses of such an approximation in the EOMCC context and discussed the lack of size intensivity of the excited state variant. It has been further shown by Peng and Kowalski [28] that in order to maintain size extensivity of the resulting Green's function, the rank of n should not exceed by more than one in comparison to m.
In this spirit, we choose the ground state cluster op-eratorT andΛ operator up to rank 2. For the (N + 1) and (N − 1) electronic problems, we choose theR and L up to excitation rank 3. This approximation not only allows us to avoid the steep n 8 b computational scaling, but also it allows for memory savings -we no longer need arrays of dimension (n 3 o n 3 v ) (n o and n v stand for number of occupied and virtual orbitals, respectively) for T 3 and Λ 3 . The maximum size of necessary arrays for R 3 and L 3 is now We can rationalize the EOM-CCSD (m,n) approximation by using a biorthogonal projection operator ( kR k eT |Φ 0 Φ 0 |L k e −T ) for the (N − 1) states in Eq. 2, resulting in the following expression: where, and ∆E k corresponds to the electron detachment energy, when R k and L k are the ket and bra eigenvectors of the N −1 electronic problem. With GFCCSDT(2,3) approximation, |Ψ N −1 k spans much larger space since it contains 3h-2p states as well. Likewise, |Ψ N +1 k also includes 3p-2h states, thus with GFCCSDT(2,3) approximation, we span larger space than GFCCSD alone. This space is necessary to describe both the ionization and attachment processes that affect the CCSD wavefunction.
The inversion present in Eq. 2 is carried out using a continued fraction formula [22]. Consequently, the computational scaling of our approach is independent of the size of the frequency grid since the Green's function at any arbitrary frequency (both on the real and imaginary frequency axis) can be evaluated using the continued fraction formula.
In order to assess the accuracy of the solvers based on GFCCSDT and GFCCSDT(2,3) developed above, we considered a cubic perovskite, SrMnO 3 , which has a G-type antiferromagnetic (AFM) ordering at lowtemperature and paramagnetic (PM) ordering at room temperature. Neél temperature for SrMnO 3 has been reported to be ∼ 233−260 K [35,36]. In this work, we focus on the high-temperature PM insulating phase, for which several photoemission experimental data [37][38][39] are available, and they predict a band gap of 1.0 to 2.3 eV.
In our previous work [13], we have shown that scGW predicts an incorrect metallic PM phase for this system. Several other more sophisticated studies, such as LDA+DMFT [40][41][42] and multitier GW +EDMFT [43], fail to open the gap in the PM phase as well. Recently, our group has reported a study of the same PM phase, where SEET(GW /ED) was successfully employed to open a gap with an outer-loop self-consistency [44]. In that study, from the evaluated orbital resolved spectra, we observed that Mn:3d and O:2p orbitals have dominant contributions to bands near the Fermi level. In the SEET(GW /ED) study, we created a series of impurity problems built from Mn:3d and O:2p orbitals. The selfenergies obtained using ED for these impurities were used to correct the GW self-energy during the self-consistency SEET cycle.
In a later study [25], we have used the same SrMnO 3 PM phase to assess the accuracy of the GFCCSD solver for the impurity problems. It revealed that for some of the impurity choices, GFCCSD predicts similarly qualitative results as ED. For example, SEET(GW /GFCCSD), employing GFCCSD as a solver for the Mn:t 2g impurity, opens a gap yielding quantitatively correct spectral functions.
However, for some other impurity choices GFCCSD was not as successful, resulting in quantitatively incorrect self-energies that affected the final spectral functions. For example, with Mn:e g and O:p π impurities, SEET with GFCCSD as a solver predicts a spurious peak near the Fermi region [25] yielding incorrect spectral functions.
In this paper, we will investigate whether inclusion of triples improves upon the GFCCSD self-energies. This improvement will be analyzed by examining self-energies of various impurities of Mn:3d and O:2p orbitals. Subsequently, we will investigate the effect of the improved impurity self-energies built from Mn:3d and O:2p orbitals as shown in Table I on the full SEET(GW /GFCCSDT) calculations.
We have assembled all the integrals for those impurity problems in an open access repository [45].
We have reported two sets of results. In the first set to enable comparison between the ED, GFCCSDT, and GFCCSDT(2,3) self-energies, we were able to carry out these calculations for A, B, D and E impurities. For the remaining impurities, F,G, and H, the ED results are impossible to obtain because of the large number of orbitals (system+bath) in such impurity problems. For these impurities the total number of orbitals exceeds 16. Note also that for all the impurities reported here, the real part of Σ is completely dynamic (frequency dependent), since we have subtracted all the static self-energy contributions.
The first set of results in which we are reporting selfenergies for impurities A, B, D, E are shown in Fig. 1. For all the impurities other than O:p σ , we observe quantitative differences when their self-energies at the GFCCSD level are compared to the ED ones. Moreover, the imaginary parts of the GFCCSD self-energies for Mn:e g and O:p π display different minima at the frequency axis in comparison to ED. After adding triples correction with our low-scaling approximated GFCCSDT(2,3) scheme, we can see that they reproduce ED self-energies with better accuracy, especially there is a noticeable improvement in the low frequency regime of ImΣ and in ReΣ. Furthermore, when employing the full GFCCSDT scheme, the self-energy for Mn:t 2g and Mn:e g impurities is further improved, as we expect. However, for O:p π impurity, Im Σ is much better reproduced around minimum with GFCCSDT(2,3) in comparison to GFCCSDT. Another crucial observation is that in both for Mn:e g and O:p π impurities, GFCCSD produces a non-causal self-energy since the GFCCSD Im Σ is positive. (A causal, physical diagonal parts of Im Σ are always negative). When triples corrections are employed, in both GFCCSDT and GFCCSDT(2,3), such a behaviour is alleviated.
The next set of impurities, C, F and G, is created by combining smaller impurities together. For each of these impurities, this results in a larger size of the combined orbital space (system+bath). Consequently, we are no longer able to compare the resulting self-energies against ED. Therefore, we consider the GFCCSDT self-energy as our reference. Here, we aim to investigate the qualitative and quantitative behaviour of the GFCCSDT(2,3) self-energy. For O:2p impurity, we observe from Fig.  2 that Re Σ behaves qualitatively similar for all the methods considered here and GFCCSDT is very similar to GFCCSDT(2,3). However, for Im Σ we observe non-causal contributions, both for GFCCSD and GFCCSDT(2,3) methods. This is reminiscent of the O:p π smaller impurity, but unlike for that impurity, here GFCCSDT(2,3) failed to yield a causal self-energy. GFCCSDT, on the other hand does not suffer from this problem. Therefore, by comparing a smaller O:p π and a combined O:p impurity, we conclude that a better performance of GFCCSDT(2,3) for O:p π was probably fortuitous (perhaps due to some error cancellation). For the other larger impurities considered, that is for Mn:3d and Mn:e g + O:p z , we observe that GFCCSDT and GFCCSDT(2,3) follow each other closely in most parts except showing some quantitative difference at low frequency regime of Im Σ. GFCCSD, on the other hand differs from other schemes containing triples, both in terms of the position of the minima and absolute values.
We analyzed the accuracy of total energy obtained from various GFCC variants while taking into account a lack of self-consistent nature of the CC theory.
Fully self-consistent Green's function methods such as Green's function self-consistent second order (GF2) or GW are approximations of the exact Luttinger-Ward functional Φ. Diagrammatically, such a functional is the sum of all bold, closed, two-particle irreducible Feynman diagrams (also called "skeleton" diagrams). In such methods ∂Φ ∂G = Σ and Green's function, G(Σ, v), is a functional of the self-energy and bare Coulomb interactions. Consequently, the self-energy and Green's function are obtained self-consistently.
The GFCC Green's function is not constructed as a self-consistent functional of G and Σ. We will analyze the error introduced due to the non-self-consistent nature of the self-energy obtained from all the GFCC variants described here. In order to do so, we will consider one traced quantity, namely the total energy for the impurity. We have evaluated the total energy in two different ways.
First, we evaluated it directly from the CC wave function as where, h q p , F q p , V pq rs are the bare one-electron, Fock and bare two-electronic Hamiltonian matrix elements in the molecular orbital (MO) basis, and γ q p and Γ rs pq are the one-body and two-body reduced density matrices respectively, obtained from CC wave functions.
Second, it was calculated directly from the Galitskii-Migdal formula [46] using Σ(ω) and G(ω) obtained from various GFCC variants considered here. In this formalism the 1-and 2-body contributions to the energy are defined as where γ(G) = −G(τ = β − ) is obtained from the correlated Green's function. For the orbital-trace we use atomic orbitals (AO) here. Since we are comparing the traced quantities, the difference in the choice of orbital basis between different approaches does not play any role.
We have summarized our results in Table II. We expect that GFCCSD and GFCCSDT(2,3) will reproduce E 1b and E 2b , individually from the CCSD wave function calculation, and that GFCCSDT will be compara- ble to CCSDT wave function results. For O:p σ impurity, both E 1b and E 2b were recovered with very good accuracy by all the methods considered. This finding can be rationalized from the fact that all these methods were very accurate for evaluating the self-energy of this impurity, as discussed in the previous section. For the O:p π impurity, we observe substantial differences between CCSD(wf) and GFCCSD; both for E 1b and E 2b . When GFCCSDT(2,3) is considered for the O:p π impurity, where 3h-2p/3p-2h projections are used for the the N-1/N+1 particle states, we see a good agreement with CCSD(wf) quantities. We believe this asserts the improved behaviour of GFCCSDT(2,3) over GFCCSD. For GFCCSDT, we hoped that the contribution of the missing projections, that is, of 4h-3p/4p-3h states will be insignificant because they typically appear as a much higher-order perturbative contribution. However, our comparison does not comply with that expectation for O:p π impurity. For Mn:e g impurity, our observations are very much the same, albeit the numerical values of differences are much smaller.
For the self-consistent SEET study with GFCC solver, we have considered MnO solid, which is a prototypical strongly correlated system with AFM ordering. We have studied this system in our previous work [25] with both the ED and GFCCSD solvers. The computational setup remained the same as in Ref. [25]. For MnO, we have observed that updating the GW self-energy via the outer loop self-consistency has very little effect unlike in the case of SrMnO 3 solid. Moreover, some of the impurities did not yield correct spectral function with GFCCSD solver in comparison to the ED solver, for example, when O:2p orbitals were considered in the impurity construction. For these two reasons, MnO will be considered as a very good test of GFCC solver performance in the rest of the discussion.
In our study of MnO with the GFCCSD solver for a case that included O:2p impurity problem, we have observed a spurious peak near the Fermi region of the spectral function (see Ref. [25]). The origin of this peak can be traced to the significant discrepancy observed in the GFCCSD self-energy when compared with the ED self- energy. This difference is prominent for the imaginary part of the self-energy.
When we carried out the full self-consistency calculation with the GFCCSDT(2,3) solver the spurious peak disappeared and the overall band gap in SEET(GW/GFCCSDT(2,3)) spectrum is very close to the spectra for SEET(GW/ED) case. However, the same self-consistent calculation  Table II. One-body (E 1b ) and two-body (E 2b ) energies in a.u. evaluated with 1-RDM and 2-RDM from different CC variants and from Galitskii-Migdal formula using G(ω) and Σ(ω) obtained from different GFCC variants.
SEET(GW/GFCCSDT) with the GFCCSDT solver failed to remove the spurious peak. We have shown both spectra from SEET(GW/GFCCSDT(2,3)) and SEET(GW/GFCCSDT) in Fig. 3. This observation is certainly counterintuitive as GFCCSDT (2,3) is only an approximation to the GFCCSDT method. In the following we will try to provide a plausible explanation of the above outcome.
We plotted the self-energy obtained from ED, GFCCSDT, GFCCSDT(2,3) and GFCCSD solvers for the O:2p impurity problem in Fig. 4. At the lowfrequency regime, both GFCCSDT and GFCCSDT (2,3) show much better agreement with the ED self-energy in comparison to the GFCCSD solver. However, if we analyze the Re Σ from all the solvers we observe that at highfrequency regime there is significant difference among the solvers -GFCCSDT(2,3) is numerically much closer to what ED predicts, and GFCCSDT is significantly different. We have plotted those quantities in the inset of Fig. 4. This large difference in static self-energy causes the chemical potential value in the first iteration of SEET with the GFCCSDT solver to be ∼79 mH lower than ED. For GFCCSDT (2,3) this difference is only ∼13 mH. We attribute this as the primary reason for the appearance of the spurious peak at the Fermi region.
We furthermore attributed the reason behind this discrepancy to the non self-consistent nature of the GFCC self-energies. To further illustrate this behavior we have evaluated a static contribution Σ ∞ to the self-energy in two different ways: 2. using the first-order self-energy diagrams. In a self-consistent approach, these two ways of evaluating static contribution to the self-energy will produce the same result (up to the numerical convergence). For example, with an ED solver, these two approaches produce the same self-energy. However, both for GFCCSDT and GFCCSDT(2,3), Σ ∞ differ quantitatively between these two ways of evaluating it. In our work so far we have taken the first approach mentioned above, where we observe a deviation from ED of ∼ 7 mH for GFCCSDT and ∼ 2 mH for GFCCSDT (2,3). When this quantity was evaluated with the second approach, surprisingly the deviation for GFCCSDT reduced to 2 mH, but for GFCCSDT(2,3) it increased to 21 mH.
We have presented a validation of the newly developed variants of the GFCC solver with triple excitations, namely GFCCSDT and GFCCSDT(2,3), on selected impurity problems created for SrMnO 3 solid. These benchmark problems demonstrated that the solver involving triples leads to both qualitative and quantitative improvement of self-energies when compared to the GFCCSD solver.
Similarly, we have analyzed the 1-and 2-body energy obtained from both CC and GFCC approaches. While for non-self-consistent approaches, such as GFCC, the energies obtained from CC alone and GFCC do not need to agree, we would expect them to be close when sufficient excitations are present in the CC wave function ansatz. Our results show that the addition of triple excitations generally improved this agreement in a substantial manner.
Finally, for MnO solid we have used the GFCCST and GFCCST(2,3) solver in combination with SEET embeding procedure on top of GW calculations. We have observed that SEET(GW/GFCCST(2,3)) gave excellent results reproducing the experimental XPS and BIS spectra. We did not achieve a similar success while employing the GFCCST solver. This unexpected result was explained by analyzing the static part of the GFCCSDT self-energy.
The inclusion of triples correction certainly improves the quality of various dynamical quantities, for example, self-energy, Green's function etc. obtained. How- ever, the non-self-consistent nature of GFCC precludes the expected systematic improvement. Perhaps a different construction of static self-energy, as alluded in the results section, will improve the final quantities obtained from a self-consistent calculation. In future, we would like to resolve this particular aspect of the theory. Moreover, we can reduce the computational complexity of the triples solver further by employing a so-called AO-driven implementation of the most expensive terms. A.S., C.-N. Y., B.P., K.K. and D. Z. acknowledge support of the Center for Scalable, Predictive methods for Excitation and Correlated phenomena (SPEC), which is funded by the U.S. Department of Energy (DOE), Office of Science, Office of Basic Energy Sciences, the Division of Chemical Sciences, Geosciences, and Biosciences.