Exploring connections between statistical mechanics and Green's functions for realistic systems. Temperature dependent electronic entropy and internal energy from a self-consistent second-order Green's function
EExploring connections between statistical mechanics and Green’s functions forrealistic systems. Temperature dependent electronic entropy and internal energy froma self-consistent second-order Green’s function.
Alicia Rae Welden, ∗ Alexander A. Rusakov, † and Dominika Zgid ‡ Department of Chemistry, University of Michigan, Ann Arbor, Michigan 48109, USA
Including finite-temperature effects from the electronic degrees of freedom in electronic structurecalculations of semiconductors and metals is desired; however, in practice it remains exceedinglydifficult when using zero-temperature methods, since these methods require an explicit evaluationof multiple excited states in order to account for any finite-temperature effects. Using a Matsub-ara Green’s function formalism remains a viable alternative, since in this formalism it is easier toinclude thermal effects and to connect the dynamic quantities such as the self-energy with staticthermodynamic quantities such as the Helmholtz energy, entropy, and internal energy. However,despite the promising properties of this formalism, little is know about the multiple solutions ofthe non-linear equations present in the self-consistent Matsubara formalism and only a few casesinvolving a full Coulomb Hamiltonian were investigated in the past. Here, to shed some light ontothe iterative nature of the Green’s function solutions, we self-consistently evaluate the thermody-namic quantities for a one-dimensional (1D) hydrogen solid at various interatomic separations andtemperatures using the self-energy approximated to second-order (GF2). At many points in thephase diagram of this system, multiple phases such as a metal and an insulator exist, and we areable to determine the most stable phase from the analysis of Helmholtz energies. Additionally, weshow the evolution of the spectrum of 1D boron nitride (BN) to demonstrate that GF2 is capableof qualitatively describing the temperature effects influencing the size of the band gap.
I. INTRODUCTION
In molecular quantum-chemical calculations of ther-modymanic properties such as Gibbs energy [1, 2], thetemperature dependent component is usually dominatedby vibrational contributions. This is due to the largegaps between electronic states, which ensure that the ex-cited state populations will be negligible. Consequently,common molecular calculations do not explicitly includetemperature effects on the electronic structure.However, for materials such as doped semiconduc-tors [3] the magnitude of the electronic band gap canbe relatively small, or for metals nonexistent altogether,allowing electronic states other than the ground stateto be accessible at low temperatures. Thus, it is nec-essary to include temperature effects into the electronicdescription. Even though for most materials the vibra-tional contribution to the specific heat is much largerthan the electronic one, there are cases when the in-corporation of the electronic contribution is necessary.The electronic contribution to specific heat is importantfor (i) heavy fermion materials at low temperatures [4], (ii) materials that do not undergo a structural transi-tion that changes the vibrational contribution, causingthe stability of phases to primarily depend on the rel-ative electronic contribution to Gibbs energy, (iii) ma-terials with a structural transition where the differencebetween vibrational contributions of the phases is of the ∗ Corresponding author: [email protected] † Corresponding author: [email protected] ‡ Corresponding author: [email protected] same order as the difference between electronic contri-butions [5, 6]. Thus, modern materials calculations canbenefit from computational tools that provide access tothe temperature dependent electronic contribution to thespecific heat, Gibbs energy, entropy, or electronic part ofthe partition function.While in traditional quantum chemical calculationsevaluating the electronic contribution to temperature de-pendent quantities is certainly not wide spread, a num-ber of such methods exist, most notably, the finite tem-perature Hartree-Fock (HF) [7], density functional the-ory (DFT) [8–11], Møller-Plesset second order perturba-tion theory (MP2) [12], coupled-cluster (CC) [13–15], andLanczos method [16, 17] for finite temperature configura-tion interaction (CI) calculations. However, these meth-ods are usually quite difficult to implement and costly touse since they rely on the modification of the parent zero-temperature method to the finite temperature formalismby adapting it to work in the canonical or grand canoni-cal ensemble. For example, to carry out such calculationsin the CI formalism, one must obtain the excited statesand corresponding Boltzmann factors in order to evalu-ate the partition function and thermodynamic averages,making application of finite temperature variants of thesemethods quite cumbersome.Conversely, for the Green’s function formalism theconnection to thermodynamics arises in a straightfor-ward manner and was derived in numerous books inthe past [18–21]. While this theoretical connection iswell understood, the actual numerical calculations of thethermodynamic quantities still remain quite challeng-ing since a fully self-consistent imaginary axis (Matsub-ara) Green’s function calculation is desired. A non-self-consistent Green’s function can result in non-unique ther- a r X i v : . [ phy s i c s . c h e m - ph ] O c t modynamic quantities.The fully self-consistent Green’s function calculationsare challenging for multiple reasons, such as large imag-inary time and frequency grids required for convergenceor multiple solutions that can be present due to the non-linear nature of the equations. It is for reasons such asthese that in recent years calculations that capitalizedon Green’s function language and yield thermodynamicquantities where mostly done for model systems [22–26].Here, we would like to stress that while multiple largescale real axis Green’s function calculations are per-formed at present for the single shot G W , GW , orsemi-self consistent GW for large realistic systems [27–33], currently, only a few research groups have managedto rigorously generalize the self-consistent finite tempera-ture (imaginary axis) Green’s function formalism to dealwith a general Hamiltonian containing all the realistic in-teractions [34–39]. Thus, any insight gained from study-ing even simple periodic systems and analyzing the possi-ble self-consistent solutions of the Matsubara formalismremains valuable.To the best of our knowledge, here, we present the firstapplication of the fully self-consistent finite temperatureGreen’s function formalism to evaluate thermodynamicquantities and phase stability for a periodic system de-scribed by a full quantum chemical Hamiltonian. Wedemonstrate that Green’s function formalism leads to asimple calculation of the electronic contribution to theHelmholtz or Gibbs energy, entropy, grand potential, andpartition function without explicitly performing any ex-cited state calculations. The presented formalism is exactat the infinite temperature limit since the perturbativeGreen’s function formulation is a perturbation that con-tains the inverse temperature as small parameter.This paper is organized as follows. In Sec. II, we intro-duce the imaginary Green’s function formalism and itsconnection to thermodynamics. In Sec. III and IV, welist properties of the Luttinger-Ward functional which isour main computational object and we explain its evalu-ation within a self-consistent Green’s function second-order (GF2) periodic implementation. In Sec. V, wepresent numerical results first for a benchmark molecularproblem and then for two 1D systems: periodic hydrogenas well as boron nitride. Finally, we form conclusions inSec. VI. II. CONNECTION WITH THERMODYNAMICS
One of the first descriptions of the connection be-tween the Green’s function formalism and thermodynam-ics was presented in the book by Abrikosov, Gorkov, andDzyaloshinski [21]. Since then multiple texts have ap-peared that discuss this connection [18, 20, 40]. For anexcellent, detailed derivation, we encourage the readerto follow Ref. 18; here, we will only mention few basicGreen’s function equations for the sake of completeness. The one-body Green’s function is defined as G ji ( z , z ) ≡ Tr[ e − β ˆ HM ˆ G ji ( z ,z )]Tr[ e − β ˆ HM ] (1)= i Tr[ τ { e − i (cid:82) γ d ¯ z ˆ H (¯ z ) ˆ d j,H ( z ) ˆ d † i,H ( z ) } ]Tr[ τ { e − i (cid:82) γ d ¯ z ˆ H (¯ z ) } ] , where ˆ d j,H ( z ) and ˆ d † i,H ( z ) are the second-quantized an-nihilation and creation operators in the Heisenberg rep-resentations and β = 1 / ( k B T ) is the inverse temperaturewhile T is the actual temperature and k B is the Boltz-mann constant. This Green’s function, depending on howthe contour γ in the complex plane is closed, can be usedto describe system’s time evolution (when z and z areset to the real-time variables), zero temperature phenom-ena, or equilibrium phenomena at finite temperature.Here, we are interested in a formalism used to calculatethe initial ensemble average that is applied to systems inthermodynamic equilibrium at finite temperature. Thisapproach is called Matsubara formalism or the “finite-temperature formalism”. For this reason, in the Green’sfunction from Eq. 1 we set z = t − iτ and z = t − iτ .Consequently, the Green’s function G Mji ( τ , τ ) = i (cid:8) θ ( τ − τ ) Tr[ e ( τ − τ − β ) ˆ HM ˆ d j e ( τ − τ
1) ˆ HM ˆ d † i ]Tr[ e − β ˆ HM ] ± θ ( τ − τ ) Tr[ e ( τ − τ − β ) ˆ HM ˆ d † i e ( τ − τ
2) ˆ HM ˆ d j ]Tr[ e − β ˆ HM ] (cid:9) (2)does not describe any time evolution of the system understudy. Instead, in this Green’s function the initial state ofthe system can be the thermodynamic state correspond-ing to a Hamiltonian, ˆ H M = ˆ H ( t ) − µ ˆ N , where ˆ N isthe particle number operator. Thus, from this Matsub-ara Green’s function, the initial ensemble average of anyone-body operator ˆ O = (cid:80) ij O ij ˆ d † i ˆ d j , can be evaluatedsimply as O = Tr[ e − β ˆ H M ˆ O ]Tr[ e − β ˆ H M ] = ± i (cid:88) ij O ij G Mji ( τ ) , (3)where τ = τ − τ since the one-body Green’s functionmatrix elements depend only on the difference of theimaginary time variables. Furthermore, the imaginarytime Green’s function G ( τ ) can be Fourier transformed tothe imaginary frequency Green’s function G ( iω n ) wherefor fermions the frequency grid is given by ω n = (2 n +1) πβ with n defined here as a positive integer. For simplicity,in the remainder of this paper, we will drop subscript M denoting Matsubara Green’s functions since from now onwe will only discuss this finite temperature formalism.The discussion above shows that at finite temperatureone can get grand canonical ensemble averages of one-body operators using the Matsubara formalism. How-ever, let us ask one more question: Can we get a sys-tem’s static thermodynamic variables such as electronicGibbs or Helmholtz energy, internal energy, and elec-tronic entropy from dynamic (frequency dependent) vari-ables such as Green’s function and self-energy?The thermodynamics of a system can be described bya thermodynamical potential, such as the grand poten-tial Ω [18]Ω = 1 β { Φ − Tr[Σ G + ln(Σ − G − )] } , (4)where the self-energy Σ = Σ( iω n ) describes all the fre-quency dependent correlational effects present in the sys-tem and G = G ( iω n ) is the reference (usually non-interacting) system’s Green’s function and G = G ( iω n )is the interacting Green’s function. The interacting andnon-interacting Green’s functions are connected throughthe Dyson equation, Σ = G − − G − . The Φ functionalfrom Eq. 4 is called the Luttinger-Ward functional[41][42] and is defined asΦ = ∞ (cid:88) m =1 m Tr[ (cid:88) n Σ ( m ) ( iω n ) G ( iω n )] , (5)where Σ ( m ) ( iω n ) is a self-energy containing all irreducibleand topologically inequivalent diagrams of order m . Thedetailed derivation of Eq. 4 is presented in Ref. 18.Thus, the computational object that provides a con-nection between the static and dynamic quantities isthe Luttinger-Ward functional. This scalar functionalΦ = ˆΦ[ G ] depends on the Green’s function and has mul-tiple important properties for Green’s function theory. III. PROPERTIES OF THE LUTTINGER-WARDFUNCTIONAL
The formal properties of Luttinger-Ward functionalhave been discussed extensively before. The functionalhas previously been applied to calculate energies foratoms and molecules [34, 35], the electron gas [43], as wellas for the Hubbard lattice [23, 24, 44, 45]. In the follow-ing section, we will only outline some of the most salientproperties of the Luttinger-Ward functional to benefitthe reader.
A. Self-energy as a functional derivative
The self-energy Σ = Σ( iω n ) can be obtained as a func-tional derivative of the Luttinger-Ward functional β δ ˆΦ[ G ] δG = ˆΣ[ G ] . (6)Here, the self-energy is defined as a functional of Green’sfunction that is evaluated independent of the Dysonequation. Consequently, in the non-interacting limit,where Σ = 0, it follows that the Luttinger-Ward func-tional is zero itself, ˆΦ[ G ] = 0. B. Connection with the grand potential Ω The grand potential is a number but the mathematicalobject defined in Eq. 4 can be viewed more generally asa functional of Green’s functions Ω[ G ]. When a Green’sfunction is a self-consistent solution of the Dyson equa-tion then the functional derivative δ Ω[ G ] δG = 0 since δ Φ =Tr[Σ δG ] and δ Ω[ G ] = β { δ Φ − Tr[ δ Σ G + Σ δG − Gδ Σ] } .Consequently, we can conclude that the functional Ω[ G ]at the stationary point is equal to the grand potential.Having grand potential one gains access to the partitionfunction (Z) since Ω = − β ln Z. (7)The Helmholtz energy, A = E − T S , where E is theinternal energy and S is the entropy of a system at agiven temperature T , is connected to the grand potentialas A = Ω + µN , where N is the number of electrons inthe system. Thus, knowing the grand potential, we caneasily calculate the electronic entropy as S = E − Ω − µNT (8)as long as we have access to the internal energy of asystem. The internal energy at a given temperature T can be evaluated using the Galitskii-Migdal formula E = 12 Tr [( h + F ) γ ]+ 2 β N ω (cid:88) n Re (cid:0) Tr[ G ( iω n )Σ( iω n )] (cid:1) , (9)where γ is the one-body density matrix, h is the one-bodyHamiltonian, F is the Fock matrix of a system, and N ω isthe size of the imaginary grid. Consequently, having ac-cess to the Luttinger-Ward functional of a system yieldsmultiple electronic thermodynamic quantities. C. Universality
Given two systems A and B at the same phys-ical temperature T and the same chemical po-tential µ described by two Hamiltonians ˆ H A = (cid:80) ij t Aij a † i a j + (cid:80) ijkl v ijkl a † i a † j a l a k and ˆ H B = (cid:80) ij t Bij a † i a j + (cid:80) ijkl v ijkl a † i a † j a l a k such that they have the same two-body integrals v ijkl but different one-body integrals t A (cid:54) = t b the Luttinger-Ward functional is same (universal) forboth of them. Since ˆΦ A = ˆΦ B , then it must also holdthat ˆΣ A ( G ) = ˆΣ B ( G ). In other words, two systems de-scribed by different G but having the same two-bodyinteractions are described by the same Luttinger-Wardfunctional. IV. EVALUATION OF LUTTINGER-WARDFUNCTIONAL WITHIN GF2A. Description of the GF2 algorithm
The self-energy, which we evaluate self-consistently inthis work, is computed perturbatively at the second-order(GF2) level. GF2 is advantageous for many reasons,namely, among others it behaves qualitatively correctfor moderately strongly correlated systems [36], unlikemethods such as MP2 or CCSD which tend to divergein these cases. GF2 has small fractional charge and frac-tional spin errors [46]. GF2 is carried out self-consistentlyon imaginary time τ and imaginary frequency iω n axeswith a computational scaling of O ( n τ N ) for molecularcases, where n τ is a prefactor that depends on the size ofthe imaginary time grid and N is the number of orbitalspresent in the problem. We build the Green’s functionusing the following expression G ( iω n ) = [( iω n + µ ) S − F − Σ( iω n )] − (10)where F and S are the Fock and overlap matrices inthe atomic orbital (AO) basis, respectively, and µ is thechemical potential, which guarantees a correct particlenumber. To obtain Σ( iω n ), we solve the Dyson equationgiven as Σ( iω n ) = G ( iω n ) − − G ( iω n ) − (11)where G ( iω n ) = [( iω n + µ ) S − F ] − is the non-interactingGreen’s function while G ( iω n ) is the interacting Green’sfunction since it contains the self-energy. To reduce thenumber of necessary grid points, we employ a spline inter-polation method to evaluate the Green’s function [38] andLegendre orthogonal polynomials to expand the Σ( τ ) ma-trix [47]. As pointed out previously, this self-consistentevaluation guarantees that both the Galitskii-Migdal andLuttinger-Ward energies are stationary with respect tothe Green’s function. For a full discussion of the algo-rithm and implementation details of GF2, we refer thereader to Refs. 36 and 46. The main computational ob-ject in our evaluation is the self-energy in the AO basis,which can be expressed asΣ ij ( τ ) = − (cid:88) klmnpq G kl ( τ ) G mn ( τ ) G pq ( τ ) ×× v ikmq (2 v ljpn − v pjln ) . (12)where v ijkl are the two-electron integrals in the AO basisin the chemist’s notation. Eq. 12 is evaluated in a Leg-endre polynomial basis [47] to accelerate the calculation.The details of a periodic GF2 implementation havebeen reported in Ref. [37]. The basic differences from themolecular version include solving the Dyson equation in k -space via G k ( iω n ) = (cid:2) ( iω n + µ ) S k − F k − Σ k ( iω n ) (cid:3) − (13) and complicating the real-space self-energy, Σ ( τ ), eval-uation by additional cell index summations:Σ ij ( τ ) = − (cid:88) g ,..., g (cid:88) klmnpq G g g k l ( τ ) G g g m n ( τ ) G g g p q ( − τ ) ×× v g g i m q k (2 v gg g g j n p l − v gg g g j l p n ) . (14)Σ k ( iω n ) in Eq. 13 and Σ ( τ ) in Eq. 14 are intercon-vertible via corresponding Fourier transforms from the k - to real-space and between the imaginary time andimaginary frequency domains also via Fourier transform.The self-energy calculation according to Eq. 14 resultsin O ( N N cell n τ ) formal scaling of the computation costwith N the number of orbitals in the unit cell and N cell the number of real space cells. This is typically a com-putational bottleneck of the GF2 self-consistency proce-dure.The level of self-consistency at which the correlatedGreen’s function equation should be iterated can dependon particular phases present in the phase diagram, suchas Mott (see Ref. [48]). In particular, the non-interactingGreen’s function G ( iω n ) build using updated Fock ma-trix or the correlated Green’s function G ( iω n ) can re-enter the evaluation of the self-energy. We observe thatthe use of G ( iω n ) with the updated Fock matrix inthe self-consistent evaluation of the self-energy is a well-behaved procedure when the strong correlations and theMott phases emerge. The full self-consistent cycle, with G ( iω n ) re-entering the evaluation of the self-energy be-comes ill behaved for these cases and we experienced diffi-culty converging it. We therefore use the former “partial”self-consistency ( G ( iω n ) with the updated Fock matrixin the self-consistent evaluation of the self-energy) for theMott regime. Such scheme is not uncommon in DMFTtype calculations [48].The expression for the total energy likewise acquires acell summation according to Eqs. 13 and 14 in Ref. [37]: E tot = E b + E b (15)= (cid:80) g ,i,j γ ij (2 h ij + [Σ ∞ ] ij ) ++ β (cid:80) g ,i,j Re (cid:104)(cid:80) n G ij ( iω n )Σ ij ( iω n ) (cid:105) . B. Evaluation of the grand potential in the k-space
The evaluation of the Luttinger-Ward functional inconserving approximations [49] such as the self-consistentsecond-order Green’s function (GF2)[36, 46, 47] is quitestraightforward and was originally derived by Luttingerand Ward [41] asΦ (2) = 14 Tr[ (cid:88) n Σ (2) ( iω n ) G ( iω n )] , (16)where Σ (2) ( iω n ) is the frequency dependent part ofthe second-order self-energy. Since both G ( iω n ) andΣ (2) ( iω n ) are readily available from a GF2 calculation,we are able to easily evaluate all terms of the functional.The expression for grand potential from Eq. 4 can beconveniently reformulated asΩ LW = 12 Tr[ γ Σ ∞ ] + Tr[ G Σ]+Tr[ln { − G Σ } ] + Tr[ln { G − } ] . (17)An excellent derivation of the above equation is givenin Ref. 35 for molecular systems. However, as we men-tioned before, for molecular systems the changes due tothe electronic contributions of the Gibbs or Helmholtzenergy are negligible due to the size of the gap.Here, we list detailed steps that need to be executedwhen dealing with crystalline systems where the elec-tronic effects influencing Helmholtz energy can be signif-icant. In the periodic implementation, we evaluate thegrand potential per unit cell, Ω . The overall expressionis given as a sum of all the componentsΩ = Ω Tr[GΣ] + Ω Tr[ln { − GΣ } ] +Ω Tr[ln { G − } ] + Ω Tr[ γ Σ ∞ ] (18)where we sum the contribution per unit cell for each termpresent in Eq. 17. Due to the crystalline symmetry, it isoften more convenient to calculate these quantities in k -space and transform the resulting quantity to the realspace rather than using an explicit real space represen-tation of all the quantities involved. All terms in Eq. 18containing the Green’s function G , or self-energy Σ werecomputed in k -space and then Fourier transformed toreal space. For example, to calculate the Ω Tr[GΣ] contri-bution, we can first evaluate the k -dependent quantityΩ k Tr[GΣ] = 2 β N ω (cid:88) i,j,n G k ( iω n ) ij Σ k ( iω n ) ji (19)where the indices i and j run over all atomic orbitalsin a given k − block. Subsequently, we perform Fouriertransform of Ω k Tr[GΣ] to yield the real space Ω Tr[GΣ] con-tribution.The most cumbersome evaluation is of the termΩ Tr[ln { − GΣ } ] , which requires diagonalization of the ma-trix G k Σ k +( G k Σ k ) † − G k Σ k ( G k Σ k ) † , where the matrices G k and Σ k are understood to be dependent on iω n . Oncethe eigenvalues of this matrix, (cid:15) ki , have been computedfor each imaginary frequency point, iω n , we haveΩ k Tr[ln { − GΣ } ] = 2 β N ω (cid:88) n,i ln { − (cid:15) ki } , (20)which can be Fourier transformed to the real space. To include the contribution from the Fock matrix,Ω Tr[ln { G − } ] , we calculateΩ k Tr[ln { G − } ] = (cid:40) β (cid:80) i ln(1 + e β ( (cid:15) ki − µ ) ) + (cid:15) ki , if (cid:15) ki − µ < β (cid:80) i ln(1 + e − β ( (cid:15) ki − µ ) ) , otherwise . (21)where by analyzing if the term (cid:15) ki − µ is smaller or greaterthan zero we account for the cases where the absolutevalue of the Fock matrix eigenvalue can be large leadingto numerical problems if only one branch of the aboveexpression is used. This term accounts for occupationchanges with temperature. For high β values (low valuesof the actual temperature T ), the expression reduces toΩ k Tr[ln { G − } ] = (cid:88) N occ (cid:15) ki , (22)which allows electrons to occupy only the lowest availablestate, as is expected at a very low temperature.We calculate the term Ω Tr[ γ Σ ∞ ] directly in the realspace as Tr[ γ Σ ∞ ], since in the AO basis, the decay ofboth γ which is the density matrix and Σ ∞ which is thefrequency independent part of the self-energy is rapidenough for a relatively few number of cells to assure aconverged value of Ω Tr[ γ Σ ∞ ] . V. RESULTSA. HF molecule
In this subsection, for a simple molecular example, ahydrogen fluoride molecule, we provide a calibration ofthe thermodynamic quantities such as internal energy E ,Helmholtz energy A , and entropy S , which are evaluatedat the GF2 level and compared to the full configurationinteraction (FCI) calculation. The evaluation of the ther-modynamic quantities at the FCI level can be done onlyfor very small molecular examples, such as HF, since sucha system has only 10 electrons and 6 basis functions inthe STO-3G basis, resulting in a small number of pos-sible configurations necessary to evaluate the FCI grandpotential. Note that to describe a true physical systemat very high temperature we would require a very largebasis set. Here, we use HF as a model molecular systemthat is calculated in a minimal basis set solely to en-able comparison of GF2 with FCI. The full configurationinteraction (FCI) quantities for HF molecule were calcu-lated previously by Kou and Hirata [50]. In addition, weshow the same system calculated with finite-temperatureat the Hartree-Fock level. We would like to emphasizethat we provide this molecular example as a benchmarkonly, in order to compare our method with highly accu-rate FCI quantum chemical data. Typically electroniccontributions to thermodynamics are not considered formolecular systems.Let us first note that in order to calculate the FCIpartition function and subsequently grand potential inthe grand canonical ensemble, we need to evaluate Z GC = n (cid:88) N =0 (cid:88) S z (cid:88) i (cid:104) Φ ( N,S z ) i | exp {− β ( ˆ H − µ ˆ N ) }| Φ ( N,S z ) i (cid:105) , (23)where the Φ ( N,S z ) i is the FCI wave function with N elec-trons and S z quantum number and the number of pos-sible occupation runs from 0 to 2 n , where n is the num-ber of orbitals. Consequently, we need to explicitly ob-tain the information about every possible excited statepresent in the system with different number of electrons.Such a task quickly becomes impossible for any largersystems. In contrast, in Green’s function methods, wenever need to explicitly evaluate any information con-cerning specific excited states. It is sufficient to evaluateEq. 17 and then Eq. 7 to obtain the grand canonical par-tition function. Thus, even for relatively large systemssuch calculations remain feasible.Results from our calibration are shown in Fig. 1. Thenumerical values used in the plots are tabularized in theSupplemental Information. The detailed description ofthe grids on which we evaluate Σ( τ ) and G ( iω n ) can befound in Ref. [51]. For the hydrogen fluoride molecule thetemperature range is huge due to the size of the Hartree-Fock HOMO-LUMO gap in this system, which is around1.0 a.u. corresponding to a temperature of around 3 . × K. Consequently, to make every state accessible tothe electrons in this system, we require extremely hightemperatures, as indicated by our results.In the very high temperature limit, both finite tem-perature Hartee-Fock and GF2 yield the internal energy(E), Helmholtz energy (A), and entropy (S) in excellentagreement with FCI results. This is of course expectedsince at very high temperatures the electronic behavior iswell described by mean field theories. For the intermedi-ate temperatures, GF2 thermodynamic quantities (E, S,and A) are closer to FCI than thermodynamic quantitiesobtained in finite temperature Hartee-Fock. For this sys-tem at intermediate temperatures, the GF2 thermody-namic quantities are always overestimated while the finitetemperature Hartree-Fock always underestimate them incomparison to FCI. The GF2 and Hartree-Fock entropiesfor low temperatures are well recovered and comparable.As expected, the low temperature GF2 internal energyis closer to FCI than the one evaluated using finite tem-perature Hartree-Fock. For our very lowest temperature(10 K), we recover a small negative entropy. This is anumerical artifact brought about from the level of con-vergence of the internal energies (1.0 × − a.u.) and theexpression for entropy which requires multiplication by β , S = β ( E − Ω − µN ). Thus, the smallest error ≈ − inthe energy will result in ≈ − error in the entropy dueto multiplication by β = 100. B. Periodic calculation of 1D hydrogen
In our previous work, GF2 was implemented for peri-odic systems and applied to a 1D hydrogen solid [37] inthe mini-Huzinaga [52] basis set. We consider the samesystem in this work, where we have used 5,000 Matsubarafrequencies to discretize the Green’s function in the fre-quency domain, 353 imaginary-time points, and 27 Leg-endre polynomials. We have found this grid size is suffi-cient to evaluate energy differences between the systemsat various temperatures. Up to 73 real space unit cellsappear in the self-energy evaluation (Eq. 14).This system is simple enough to be a test bed for self-consistent Green’s function theory; however, it displaysa phase diagram that is characteristic of realistic solids.At different internuclear separations, corresponding todifferent pressures, we were able to recover multiple so-lutions. Although yielding different electronic energiesand different spectra, these solutions can be mathemat-ical artifacts of the nonlinear self-consistency procedurepresent in GF2; however, they can also have physicalmeaning corresponding to different solid phases.To decide which phase is more stable at a given tem-perature, it is necessary to consider the Helmholtz en-ergies that we are able to obtain from the Luttinger-Ward functional for every solution. Previously, for theinverse temperature of β = 100, at most of the geome-try points, we have identified two possible phases withdifferent internal energies, E , that were obtained start-ing the iterative GF2 procedure either from an insulat-ing or metallic solution. The results of our investigationcan be found in Fig. 4 of Ref. 37. Currently, to ana-lyze the stability of the solutions, for a range of inversetemperatures β = 25 , , E , Helmholtz energy A , and the entropic contribution T S to the Helmholtz energy. We also improved our con-vergence criteria not only converging the internal energy E = E b + E b (as we have done in the previous work)but also converging both the E b , E b , and the Helmholtzenergy separately. This much more stringent procedureto analyze convergence of GF2 leads us to slightly revisedsolutions for the 1D hydrogen solid which we discuss inthe subsequent sections.The spectra are produced from analytical continua-tion of the imaginary axis Green’s function G ( k , iω n ) tothe real axis G ( k , ω n ) [53]. The spectral weight is pro-portional to Im G ( k , ω n ). A 2D color projection of thespectral function on the ( k , ω n ) plane can be viewed asa “correlated band structure” analogous to the conven-tional band structure within effective one-electron mod-els such as Hartree–Fock and DFT. As in one-electronmodels, zero spectral weight at the Fermi energy ω F isindicative of a gapped system. The peaks emerging im-mediately below and above ω F for a given k correspondto the energies of the highest occupied (HOCO) and low-est unoccupied (LUCO) crystalline orbitals, respectively.We should stress, however, that since G ( k , ω n ) is a many-body correlated Green’s function, such correspondence is T (K) − . − . − . . . . E F C I − E m e t h o d ( a . u . ) GF2HF T (K) − . − . − . − . − . . . . . A F C I − A m e t h o d ( a . u . ) GF2HF T (K) − . − . . . . . . . . S F C I − S m e t h o d ( k B ) GF2HF
FIG. 1. The differences between GF2, finite temperature Hartree-Fock and FCI for the hydrogen fluoride molecule at varioustemperatures. not rigorous and merely serves as a convenient analogy.
1. Short bond length/high pressure
At the interatomic separation of 0.75 ˚A, we recoveredonly one gapless, metallic solution for all the values of in-verse temperature ( β = 25, 75, 100). We established thatstarting from two different initial guesses leads in bothcases to two final solutions that were different in internalenergy and Helmoltz energy by less than 10 − a.u. Con-sequently, we deemed that we obtained the same metal-lic solution in both cases. The spectral functions andspectral projections for this metallic solution at differentvalues of inverse temperature are shown in Fig. 2. Forall the temperatures examined in this short bond lengthregime, the self-energy displays a Fermi liquid character.In our previous work [37], we observed a small internalenergy difference between the solutions obtained usingdifferent starting point at β = 100. We currently observethat this difference can be eliminated if we assure thatthe convergence criteria are not only fulfilled for the totalenergy E tot = E b + E b but also both the E b and E b components separately.
2. Intermediate bond length/intermediate pressure
Spectral functions and projections for the 1D hydrogensolid with an interatomic separation of 1.75 ˚A are pre-sented in Fig. 3. We have displayed differences in internalenergy (E), entropy (written as -T∆S), and Helmholtzenergy (A) in Table I. For this system at the inversetemperatures of β = 100 and 75 we obtained two solu-tions from two different initial guesses. We are able tocharacterize the first solution (“solution 1”) as a bandinsulator since the spectral function shows a gap and theself-energy displays a Fermi liquid profile. The secondsolution (“solution 2”) is gapless and therefore a metal.At an inverse temperature of β = 25, both “solution1” and “solution 2” obtained from two different startingguesses have the same spectra and identical Helmholtz β ∆ E -T ∆ S ∆ A100 -0.00370 0.10854 0.1048475 -0.00528 0.09971 0.0944325 1.36 × − -2.94 × − -2.93 × − TABLE I. Thermodynamic data for a 1D periodic hydrogensolid with separation R=1.75 ˚A. The units for β are 1/a.u.The units for all other quantities are a.u. All values are ob-tained by subtracting “solution 1” from “solution 2”. ∆ E=E sol -E sol , ∆ A =A sol -A sol , ∆ S =S sol -S sol . Note thatthe quantities at β =25 are below the precision of convergence(1 × − ). energy, indicating that there is only one stable solution— a single phase. Note that in Fig. 3 for β = 25 both thespectral functions for “solution 1” and “solution 2” seemto have different heights; however, it is an illusion sinceboth spectral functions are plotted with a different z-axisrange. We deem that the Helmholtz energy is identicalfor both these solutions at β = 25 since the obtaineddifferences are below our convergence threshold which is1 × − a.u.For two lower temperatures ( β = 100 and 75) by com-paring thermodynamic quantities, we are able to deter-mine which phase, “solution 1” or “solution 2”, is themost thermodynamically stable. From the data in Ta-ble I, we are able to determine that “solution 1” is themost stable phase from the positive value of ∆A at both β =100 and β =75. It is also interesting to note, that thissolution is stable due to the entropic factor not due tothe difference in the internal energy.
3. Long bond length/low pressure
Similar to the R=1.75 ˚A regime, for an interatomicseparation of 2.5 ˚A at inverse temperatures β = 100 and75, we recover two solutions from the two different ini-tial guesses, see Fig. 4. At a high temperature β =25 sln 1 - 0.75 A = 100 = 75 = 25 FIG. 2. Spectral functions and projections at different inverse temperatures for the metallic solution of 1D periodic hydrogensolid at R=0.75 ˚A. we observe only one solution independent of the initialguess, and thus only a single phase is present. From thepositive value of ∆A we are able to see that “solution1” is the most stable phase at β =75 and β =100 (Ta-ble II). This solution is a band insulator since it hasa Fermi liquid self-energy profile at these temperatures.The other solution, denoted as “solution 2” at low tem-perature is metallic and changes into a Mott insulator athigh temperatures. The internuclear distance of R=2.5˚A is close to the region where the phase transition occurs,thus results obtained from GF2 which is a low order per-turbation expansion may not be reliable. The low levelperturbation theories such as GF2 are known to be moreaccurate deep within a phase and can experience prob-lems near the phase transition point.Finally, at an interatomic separation of R=4.0 ˚A, forall the temperatures, both initial guesses yield the sameconverged GF2 result - a single phase. The spectra asa function of inverse temperature can be seen in Fig. 5.This single solution is a Mott insulator as confirmed bythe divergent imaginary part of the self-energy.
4. GF2 phase diagram for 1D hydrogen solid
It is instructive now to collect all the data and con-struct a simple phase diagram for the 1D hydrogen solidas a function of inverse temperature β and intermoleculardistance R. We have plotted the phase diagram in Fig. 6. β ∆ E -T ∆ S ∆ A100 0.08489 0.16257 0.2474675 0.07103 0.16155 0.2325925 6.09 × − × − × − TABLE II. Thermodynamic data for a 1D periodic hydrogensolid with separation R=2.5 ˚A. The units for β are 1/a.u. Theunits for all other quantities are a.u. All values are obtainedby subtracting “solution 1” from “solution 2”. ∆ E =E sol -E sol , ∆ A =A sol -A sol , ∆ S =S sol -S sol . Note that thequantities at β =25 are below or comparable with the precisionof convergence (1 × − ). On this diagram, for these regions where two phases co-exist, we denoted the most stable phase according to theHelmholtz energy A by framing its symbol using a blackline.For the shortest bond length (R=0.75 ˚A), the 1D hy-drogen solid remains metallic at all temperatures con-sidered and only one phase is present. At intermediatebond lengths (R=1.75, 2.0 ˚A), multiple phases coexist.At lower temperatures, we recover both a metal and bandinsulator as possible solutions, with the band insulatorbeing the most stable phase according to the Helmholtzenergies. This is in line with physical intuition that weshould recover an insulator rather than a metal at low = 100 = 75 = 25 = 25 = 75 = 100 sln2-1.75 FIG. 3. Spectral functions and projections at different inverse temperatures for the band insulator (“solution 1”) and themetallic solution (“solution 2”) of 1D periodic hydrogen solid at R=1.75 ˚A. Note the difference in scale between “solution 1”and “solution 2”. temperature. For R=1.75 ˚A, at temperatures higher than β = 50, we recover only the metallic phase.For R=2.0 ˚A at high temperature ( β = 25), we seeboth a metal and a Mott insulator coexisting, with the Mott insulator being the most stable phase. At lowertemperatures, we see the coexistance of a metallic andband insulator solution. The Helmholtz energy indicatesthat band insulator is the most stable phase.0 sln 1 - 2.5 = 100 = 75 = 25 sln2 - 2.5-correct = 100 = 75 = 25 FIG. 4. Spectral functions and projections at different inverse temperatures for the band insulator (“solution 1”) and themetallic solution (“solution 2”) of 1D periodic hydrogen solid at R=2.5 ˚A.
At a longer bond length of (R=2.5 ˚A), multiple phasesare still present. As in the cases of intermediate bondlength, we recover both a band insulator and a metal so-lution, with Helmholtz energy favoring the band solution.We would like to reiterate that a phase transition occurs somewhere in this intermediate region, and it is likelythat the results of a second-order perturbation theorymay not be accurate enough. At higher temperatures,we recover a Mott insulator as the only phase present.For the largest separation (R=4.0 ˚A), the system re-1 = 100 = 75 = 25 FIG. 5. Spectral functions and projections at different inverse temperatures for the Mott insulator present in 1D periodichydrogen solid at R=4.0 ˚A. R FIG. 6. A phase diagram containing all distances and tem-peratures calculated for a 1D hydrogen solid. Where multiplephases exist, the most stable phase is outlined in black. β isinverse temperature in units 1/a.u. and R is the separationbetween hydrogens in ˚A. mains a Mott insulator at all studied temperatures. C. Periodic calculation of 1D boron nitride
In this section we present a periodic calculation of 1Dboron nitride (BN) at inverse temperatures β =70, 75,and 100. The B–N distance is set to 1.445 ˚A, the bondlength in the corresponding 2D system [54]. For thesecalculations we used a modification of the ANO-pVDZGaussian basis set [55] from which, in order to avoidlinear dependencies, we removed the diffuse (below 0.1)exponents and polarization functions. Larger grids of30000 Matsubara frequencies and 100 Legendre polyno-mials were found necessary for an adequate discretizationfor the self-energy and Green’s function. The self-energyis evaluated encompassing 39 cells — as many as neededfor a converged Hartree–Fock exchange.Shown in Fig. 7 is the spectral projection of the real-frequency correlated Green’s function at β =100. Usingthe analogy previously discussed at the end of SectionV. B, we plot the “highest occupied” and “lowest unoc-cupied correlated bands” of 1D BN separated by a sizablegap of approximately 3.5 eV at k=- π . This magnitude ofthe band gap is expected based on the properties of its2D counterpart [54].Once again, we treat such a simple system as a bench-mark example and we are interested in the evolution ofthe 1D BN band gap as the temperature changes andinfluences electronic degrees of freedom. Let us stressthat while an evolution of the spectrum as a function oftemperature is expected, in order to reproduce it reli-2 = 100 FIG. 7. Spectral projection in the vicinity of the Fermi energyat β =100 of periodic 1D boron nitride solid. Only the “corre-lated bands” corresponding to the conventional HOCO’s andLUCO’s (see Sec. V. B) are displayed. Note that the energyscale is in eV.FIG. 8. Periodic 1D boron nitride solid spectrum for the“correlated HOCO and LUCO” (see Sec. V. B) for severaltemperatures at k= − π . The chemical potential was adjustedfor ω = 0 to fall in the middle of the band gap. ably, the computational procedure must be very robust.We not only must be able to iteratively converge Green’sfunction and self-energy yielding different Green’s func-tions at different temperatures, we also must continuethe results obtained on the imaginary axis to the realaxis risking that the continuation will obscure the spec-tral features and the gap temperature dependence willno longer be visible.In Fig. 8 we plot the spectrum at k=- π for β =100, 75,and 70. As expected, with the increase of temperature(decrease of β ) the spectral peaks corresponding to thebands broaden and the band gap reduces. Note that bothFig. 7 and 8 use eV as energy units for clarity. VI. CONCLUSION
While the theory describing the connection betweenthe Matsubara Green’s function formalism and thermo-dynamics has been known since the 1960s, few compu-tational methods are currently capable of employing theMatsubara formalism for realistic systems. This is dueto its computationally demanding nature that requirescomplicated time and frequency grids as well as the it-erative nature of the equations. Moreover, little is knowabout obtaining different solutions corresponding to dif-ferent phases or non-physical solutions that can appear asthe result of the non-linear procedure used when Green’sfunctions are constructed iteratively on the imaginaryaxis.In our last work [37], we have demonstrated one ofthe first applications of the fully self-consistent Matsub-ara Green’s function formalism to a benchmark peri-odic problem with realistic interactions, that is, a 1Dhydrogen solid. In the current work, we have shownthat it is possible to evaluate temperature dependentthermodynamic quantities using a self-consistent second-order Green’s function method. Using the self-consistentGreen’s function and self-energy, we were able to eval-uate the Luttinger-Ward functional at various tempera-tures and obtain static thermodynamic quantities suchas Helmholtz energy, internal energy, and entropy. Eval-uation of these quantities gives us access to the parti-tion function and any thermodynamic quantity that canbe derived from it. Unlike the FCI case, in GF2 wedo not need to explicitly calculate excited state energiesor Boltzmann factors, making calculation of thermody-namic quantities for larger systems feasible.To calibrate the thermodynamic data obtained fromGF2 against FCI, we have illustrated that for a hydro-gen fluoride molecule at high temperature, we are able toobtain energies in excellent agreement with finite temper-ature FCI. In the lower temperature regime, we observed,as expected, some deviation from the FCI answer but theoverall quality of the results still remained high.Finally, since the thermodynamic data can be used toconstruct phase diagrams, we used a simple 1D hydrogensolid to investigate possible phases at different temper-atures and interatomic distances. We obtained differentphases such as a band insulator, Mott insulator, or metal,and we were able to distinguish which phase is more sta-ble at various temperatures. Determining the stability ispossible due to our ability to calculate not only the in-ternal energy but also the entropic contribution at eachtemperature. Consequently, based on the difference ofthe Helmholtz energy, we are able to determine the moststable phase for each of the temperature points.Additionally, we have performed a calculation of 1Dperiodic BN solid demonstrating that GF2 can reproducethe gap deviation as a function of temperature. Thus, forhigher temperatures, we have observed the narrowing ofthe electronic bandgap.While we acknowledge that GF2 is a low order pertur-3bative method and as such can deliver results that areinaccurate, we believe that for most weakly and mod-erately correlated systems such as semiconductors withsmall band gaps, it can be used in the future to provideaccurate thermodynamic data and phase diagrams. Forsystems where the correlations are strong, GF2 can becombined with methods such as self-energy embeddingtheory (SEET) [39, 56, 57] to provide accurate answersconcerning thermodynamics.
ACKNOWLEDGMENTS
D.Z. and A.R.W. would like to acknowledge a NationalScience Foundation (NSF) grant No. CHE-1453894.A.A.R. acknowledges a Department of Energy (DOE)grant No. ER16391 and Prof. Emanuel Gull for help-ful discussions.
VII. APPENDIX: HIGH-FREQUENCYEXPANSION OF THE GREEN’S FUNCTIONFOR EVALUATING THE LUTTINGER-WARDFUNCTIONAL
Here, we consider the contribution to the Luttinger-Ward functional in the high frequency limit. In the Mat-subara formalism for large frequencies, the Green’s func-tion and self-energy can be expressed as a series G ( iω n ) = G iω n + G ( iω n ) + O ( 1( iω n ) ) , (24)Σ( iω n ) = Σ iω n + Σ ( iω n ) + O ( 1( iω n ) ) . (25)The coefficients for the Green’s function expansionfor non-orthogonal orbitals with a quantum chemistryHamiltonian [58] are given as G = S − ,G = S − ( F − µS ) S − . (26)The Σ and Σ coefficient of the self-energy has a com-plicated explicit form as demonstrated in Ref. 58 and isevaluated numerically asΣ = Re (Σ( iω max ) × iω max )Σ = (cid:18) Σ( iω max ) − Σ iω max (cid:19) × ( iω max ) . (27)For the molecular example used in this work, the highfrequency contribution to Tr[ G Σ] was evaluated in thesame manner as previously discussed in Ref. 36 and 38.At the very high temperatures considered for molecular examples, the high frequency contribution becomes neg-ligible. However, it is necessary to include at the lowertemperatures considered. The high-frequency contribu-tion to the Ω kT r [ G Σ] term of Eq. 19 can be evaluated in away analogous to the molecular case.To evaluate the high frequency contribution to the log-arithm term Tr[ln { − G Σ } ] in Eq. 17, we expanded thelogarithm as a Taylor series ln(1 − x ) = x − x + . . . toyield ln (cid:16) − G Σ ( iω n ) (cid:17) = G Σ ( iω n ) + O (cid:16) iω n ) (cid:17) . (28)We exclude in this expansion all terms that are equal orsmaller in magnitude than O (cid:16) iω n ) (cid:17) . Thus, it is onlynecessary practically to evaluate the first term in the ex-pansion to capture the important contribution from thehigh frequency limit. We would like to emphasize thatalthough we do not present results including the high fre-quency contribution to the logarithm term Tr[ln { − G Σ } ]in Eq. 17 for the periodic hydrogen solid, this term wasevaluated and was found to be minuscule compared tothe magnitude of the other energies. However, we wouldlike stress that it is possible that this contribution maybe substantial for other systems. VIII. SUPPLEMENTAL INFORMATION k B )T (K) Hartree-Fock GF2 FCI Hartree-Fock GF2 FCI Hartree-Fock GF2 FCI10 -98.571 -98.588 -98.597 -98.571 -98.585 -99.845 0.0 -0.003 0.010 -98.571 -98.588 -98.597 -98.571 -98.621 -99.944 0.0 0.0 0.010 -97.944 -98.135 -98.049 -101.021 -103.067 -102.107 3.175 3.566 3.47510 -96.794 -96.988 -96.945 -150.563 -151.410 -151.244 4.979 4.949 4.95810 -92.028 -92.057 -92.056 -729.937 -730.100 -730.095 5.348 5.348 5.34810 -88.483 -88.487 -88.487 -6846.975 -6847.001 -6847.003 5.406 5.406 5.406 TABLE III. FCI, GF2, and finite temperature Hartree-Fock results for the hydrogen fluoride molecule in STO-3G basis. Ahigher value of temperature corresponds to a smaller value of β = 1 /k B T . Note that the grand potential depends on chemicalpotential Ω = E − T S − µN that for gaped system is not uniquely defined. The difference in the FCI, finite temperatureHartree-Fock, and GF2 for the 10 − K temperatures are the result of this non-uniqueness of the chemical potential. TheHelmholtz energy A = Ω + µN = E − T S , which is the sum of the grand potential and µN term does not suffer from thisnon-uniqueness. [1] D. A. McQuarrie and J. D. Simon, Molecular thermody-namics (University Science Books Sausalito, CA, 1999).[2] J. W. Ochterski, Gaussian Inc, Pittsburgh, PA , 1 (2000).[3] J. Y. Tsao,
The World of Compound Semiconductors (Citeseer).[4] E. G. Moroni, G. Grimvall, and T. Jarlborg, Physicalreview letters , 2758 (1996).[5] T.-Q. Yu and M. E. Tuckerman, Physical review letters , 015701 (2011).[6] Q. Zhu, A. G. Shtukenberg, D. J. Carter, T.-Q. Yu,J. Yang, M. Chen, P. Raiteri, A. R. Oganov, B. Pokroy,I. Polishchuk, et al. , Journal of the American ChemicalSociety , 4881 (2016).[7] N. D. Mermin, Annals of Physics , 99 (1963).[8] N. D. Mermin, Physical Review , A1441 (1965).[9] M. Dharma-wardana, arXiv preprint arXiv:1602.04734(2016).[10] J. Smith, A. Pribram-Jones, and K. Burke, arXivpreprint arXiv:1509.03097 (2015).[11] V. V. Karasiev, T. Sjostrom, J. Dufty, and S. Trickey,Physical review letters , 076403 (2014).[12] S. Hirata and X. He, The Journal of chemical physics , 204112 (2013).[13] M. Altenbokum, K. Emrich, H. K¨ummel, andJ. Zabolitzky, in Condensed matter theories (Springer,1987) pp. 389–396.[14] S. H. Mandal, R. Ghosh, G. Sanyal, and D. Mukher-jee, International Journal of Modern Physics B , 5367(2003).[15] M. R. Hermes and S. Hirata, The Journal of ChemicalPhysics , 102818 (2015).[16] J. Jakliˇc and P. Prelovˇsek, Physical Review B , 5065(1994).[17] J. Prelovsek, P. Bonca, arXiv preprint arXiv:1111.5931(2011).[18] G. Stefanucci and R. van Leeuwen, NonequilibriumMany-Body Theory of Quantum Systems: A Modern In-troduction (Cambridge University Press, 2013).[19] R. D. Mattuck,
A Guide to Feynman Diagrams in theMany-Body Problem (Courier Corporation, 2012).[20] A. L. Fetter and J. D. Walecka,
Quantum theory of many-particle systems (Courier Corporation, 2003).[21] I. Dzyaloshinski et al. , Methods of quantum field theoryin statistical physics (Courier Corporation, 1975).[22] G. Li, W. Hanke, A. N. Rubtsov, S. B¨ase, and M. Pot-thoff, Physical Review B , 195118 (2009).[23] M. Potthoff, M. Aichhorn, and C. Dahnken, PhysicalReview Letters , 206402 (2003).[24] M. Potthoff, The European Physical Journal B-Condensed Matter and Complex Systems , 335 (2003).[25] M. Potthoff, arXiv:1407.4065 [cond-mat.str-el].[26] M. Potthoff, The European Physical Journal B-Condensed Matter and Complex Systems , 429 (2003).[27] D. Neuhauser, Y. Gao, C. Arntsen, C. Karshenas, E. Ra-bani, and R. Baer, Physical review letters , 076402(2014).[28] L. Hung, H. Felipe, J. Souto-Casares, J. R. Chelikowsky,S. G. Louie, and S. ¨O˘g¨ut, Physical Review B , 085125(2016).[29] M. van Schilfgaarde, T. Kotani, and S. Faleev, Physicalreview letters , 226402 (2006). [30] F. Caruso, P. Rinke, X. Ren, A. Rubio, and M. Scheffler,Physical Review B , 075105 (2013).[31] S. V. Faleev, M. Van Schilfgaarde, and T. Kotani, Phys-ical review letters , 126406 (2004).[32] C. Rostgaard, K. W. Jacobsen, and K. S. Thygesen,Physical Review B , 085103 (2010).[33] M. Govoni and G. Galli, Journal of chemical theory andcomputation , 2680 (2015).[34] R. van Leeuwen, N. E. Dahlen, and A. Stan, PhysicalReview B , 195105 (2006).[35] N. E. Dahlen, R. van Leeuwen, and U. von Barth, Phys-ical Review A , 012511 (2006).[36] J. J. Phillips and D. Zgid, The Journal of ChemicalPhysics , 241101 (2014).[37] A. A. Rusakov and D. Zgid, The Journal of chemicalphysics , 054106 (2016).[38] A. A. Kananenka, A. R. Welden, T. N. Lan, E. Gull, andD. Zgid, Journal of Chemical Theory and Computation(2016).[39] A. A. Kananenka, E. Gull, and D. Zgid, Physical ReviewB , 121111 (2015).[40] W. H. Dickhoff and D. Van Neck, Many-body theory ex-posed!: propagator description of quantum mechanics inmany-body systems (World Scientific, 2008).[41] J. M. Luttinger and J. C. Ward, Physical Review ,1417 (1960).[42] (), in some texts the functional Φ is simply referred toas the Phi functional and the functional Ω is called theLuttinger-Ward functional.[43] C.-O. Almbladh, U. V. Barth, and R. v. Leeuwen, Inter-national Journal of Modern Physics B , 535 (1999).[44] V. Janiˇs, Phys. Rev. B , 11345 (1999).[45] V. Janis, Journal of Physics: Condensed Matter , L311(2003).[46] J. J. Phillips, A. A. Kananenka, and D. Zgid, The Jour-nal of chemical physics , 194108 (2015).[47] A. A. Kananenka, J. J. Phillips, and D. Zgid, Journal ofchemical theory and computation (2016).[48] S. Choi, A. Kutepov, K. Haule, M. van Schilfgaarde, andG. Kotliar, NPJ Quantum Materials , 16001 (2016).[49] G. Baym and L. P. Kadanoff, Physical Review , 287(1961).[50] Z. Kou and S. Hirata, Theoretical Chemistry Accounts , 1 (2014).[51] (), in our molecular calculations we have used 200 Leg-endre polynomials, 50,000 points on the imaginary fre-quency grid, and 5,200 imaginary-time points, corre-sponding to the full grid size before the spline interpola-tion technique was employed for the imaginary frequencygrid. One can expect a reduction to 10% of the originalgrid size using this technique. We used the most strictcriteria for the accuracy, corresponding to a δ = 10 − .Additionally, the Legendre basis vastly reduces the costarising from the number of τ points on the time grid.[52] S. Huzinaga, Gaussian basis sets for molecular calcula-tions (Elsevier, Amsterdam, 1984).[53] B. Bauer, L. D. Carr, H. G. Evertz, A. Feiguin,J. Freire, S. Fuchs, L. Gamper, J. Gukelberger, E. Gull,S. Guertler, A. Hehn, R. Igarashi, S. V. Isakov, D. Koop,P. N. Ma, P. Mates, H. Matsuo, O. Parcollet, G. Pa-wowski, J. D. Picon, L. Pollet, E. Santos, V. W. Scarola, U. Schollwck, C. Silva, B. Surer, S. Todo, S. Trebst,M. Troyer, M. L. Wall, P. Werner, and S. Wessel, Journalof Statistical Mechanics: Theory and Experiment ,P05001 (2011).[54] P. Miro, M. Audiffred, and T. Heine, Chem. Soc. Rev. , 6537 (2014).[55] F. Neese and E. F. Valeev, Journal of Chemical The-ory and Computation , 33 (2011), pMID: 26606216, http://dx.doi.org/10.1021/ct100396y.[56] T. N. Lan, A. A. Kananenka, and D. Zgid, The Journalof chemical physics , 241102 (2015).[57] T. Nguyen Lan, A. A. Kananenka, and D. Zgid, Journalof Chemical Theory and Computation , 4856 (2016).[58] A. A. Rusakov, J. J. Phillips, and D. Zgid, The Journalof chemical physics141