Power laws used to extrapolate the coupled cluster correlation energy to the thermodynamic limit
AAdvanced materials modelling in the thermodynamic limit: An analysis of theasymptotic behavior of finite size errors using coupled cluster theory
Tina N Mihm ( a ) , ( b ) , Bingdi Yang ( a ) , and James J. Shepherd ( a ) , ( b ) ∗ ( a ) Department of Chemistry, University of Iowa ( b ) University of Iowa Informatics Initiative, University of Iowa (Dated: July 24, 2020)We show that the conventional power law of the inverse of the particle number ( N − ) for ex-trapolation of total energies into the thermodynamic limit does not reproduce the Gell-Mann andBrueckner results for the three-dimensional high density electron gas, even when system sizes of N ∼ N − / fits the data beyond 1000 electrons,which has also been found by other authors. The N − / power law also appears to fit data fromother densities better than N − . To explain this observation, we develop a novel analysis of thelow- G components of the wavefunction, and show they are not well-sampled enough, even at highelectron numbers, to yield a straight-forward N − convergence to the thermodynamic limit (TDL).We show that, for most practical calculations, the TDL correlation energy lies between the extrapo-lation from N − and the extrapolation from N − / . Furthermore, the energy extrapolation from an N − / power law tends to either increase or decrease with increasing system size, but when energiesincrease is an indication that N − can effectively be used for extrapolation. Introduction:-
Coupled cluster theory, with its bal-ance of cost and accuracy in treating electronic struc-ture, has recently been highlighted for its ability to ac-curately determine defect energies [1], calculate energydifferences between spin states [2] and phases [3], iden-tify plasmons [4], calculate optical gaps [5], and describethe interactions between water and graphene [6]. Each ofthese has required a degree of detail in the wavefunctionor a correlation energy calculation ( E corr = E total − E HF ).Based on these initial applications, coupled cluster couldtransform materials design. Since coupled cluster calcu-lations are normally performed in a supercell, energiesand properties need to be carefully extrapolated to thethermodynamic, or infinite-size, limit [7–10]. There is,therefore, a critical need to understand the exact way inwhich the correlation energy approaches the thermody-namic limit, in terms of a limiting power law. As thereis some significant disagreement in the modern literatureabout what form the power law ought to take, [7–10]practitioners can be uncertain over the accuracy of ex-trapolated energies. This is important to resolve, as cou-pled cluster theory should reproduce appropriate conver-gence to the thermodynamic limit for a range of simplesystems. The rationale for this is that the random phaseapproximation (RPA), an infinite-order resummation ofring diagrams, is sufficient to reproduce the long-rangepair correlation function [10, 11] and that coupled clus-ter theory contains the appropriate diagrams to do thesame [12–15].In order to make progress towards resolving this de-bate, we set out to find the exact asymptotic power lawobeyed by the correlation energy. Our initial analysistakes advantage of the fact that the random phase ap- ∗ [email protected] proximation (RPA) is exact in the high-density limit ofthe uniform electron gas [16, 17] and that, therefore,finite-sized coupled cluster doubles (CCD) calculationswill be exact at small N [13, 14, 18]. We thereforeperformed CCD calculations of up to N = 2042 withtwist averaging. We found that an N − / , instead ofa N − , power law reproduces the exact thermodynamiclimit. This is significant, because there are a growingnumber of studies that have generally taken N − as thecorrect power law for the coupled cluster correlation en-ergy [2, 6, 14, 19] (inspired by a similar extrapolation forthe total energy).Since our findings were unexpected, we turned to ana-lytically modelling the long-range part of the wavefunc-tion, which gives rise to finite-size effects. We found thatthe there is a cross-over between the effectiveness usingan N − / power law extrapolation and an N − power lawextrapolation: the way in which the extrapolated ener-gies from the N − / power law behave with increasingsystem size determines which power law is most effec-tive. We here compare our findings with recent reportsof various power laws in the approximate ab initio cal-culations of solids [7, 19–22] and discuss the implicationsof our work for practice. Coupled cluster doubles on the uniform electrongas:-
We follow methods and procedures outlined morethoroughly elsewhere [14, 18, 23]; however, for the benefitof the reader, we will highlight some key methodologicaldetails here. Coupled cluster theory uses an exponentialansatz for the wavefunction: Ψ = e T Ψ HF .This formulation yields self-consistent equations forthe individual elements t ijab of the operator T , and anexpression for the correlation energy E corr ≡ E total − E HF as follows: E corr = 14 (cid:88) ijab t ijab ¯ v ijab (1) a r X i v : . [ phy s i c s . c o m p - ph ] J u l where ¯ v ijab are antisymmetrized four-index electron re-pulsion integrals. Here, i and j refer to occupied orbitals,and a and b refer to virtual orbitals over some finite basisset. This energy expression is related to the transitionstructure factor recently coined by Gr¨uneis and cowork-ers, S G = (cid:80) (cid:48) ijab (2 t ijab − t jiab ), where the prime indicatesthat the sum is only taken over excitations that are re-lated by the momentum transfer G [6, 19, 24].In the electron gas, all orbitals are plane waves: φ j ∝ exp( i k j · r ), where k j is a vector of the momentum quan-tum numbers in 3D, r is the electron coordinate, and i = √−
1. A simple cubic three-dimensional box is usedwith electron numbers which correspond to closed-shellconfigurations at the Γ-point. We use the Ewald in-teraction (per convention) for calculations with periodicboundary conditions. This gives rise to matrix elements v ijab (electron repulsion integrals) that take the familiarform 1 /G , where G is the magnitude of the momen-tum transfer during the i , j to a , b excitation, providedthe excitation is allowed by momentum conservation (i.e., k a − k i = k j − k b = G ). We explicitly calculate and in-clude a Madelung term, v M in our calculations [25]. Forexample, at a Wigner-Seitz radius r s = 1 . N = 14, v M = 0 . M orbitals that lie inside akinetic energy cutoff E cut,m = k .The Hartree–Fock eigenvalues for the occupied and vir-tual orbitals are, following our previous work [23]: (cid:15) p = (cid:40) k p − (cid:80) j ∈ occ p (cid:54) = j v pj − v M , p ∈ occ k p − (cid:80) j ∈ occ v pj , p ∈ virt . (2)In these equations k p is the kinetic energy of an elec-tron in orbital p and v pj is the exchange energy betweenorbitals p and j . We note that the occupied orbitals areeach lowered in energy by v M , the Madelung constant.This represents the exchange interaction of an electronwith itself. In the thermodynamic limit, v M → N and thebasis set size M ) that allows for a single CCD calcula-tion to closely approximate (within 1 mHa/el) the twist-averaged CCD energy—reducing the number of CCD cal-culations that need to be run by a factor of 100 [27]. Thesecond is the ability to extrapolate to the complete basisset and thermodynamic limits independently from oneanother [28], which we further adapt here. Connectivity twist averaging:-
Twist averaging isemployed to reduce the fluctuations in the energy withrespect to system size on converging to the thermody-namic limit [6, 9, 19, 29–35]. In general, this amounts tooffsetting the grid of momenta by applying a twist angle, k s , to each orbital such that: φ j ∝ exp( i ( k j − k s ) · r ), and then averaging the correlation energy over all twistsangles such that: (cid:104) E corr (cid:105) k s = 1 N s N s (cid:88) t =1 E corr ( k s t ) (3)where N s is the number of twist angles used in the twistaveraging. Instead of calculating this costly sum ex-plicitly, we instead use the connectivity twist averag-ing scheme. This method estimates a single twist an-gle for each calculation that has the same energy as thetwist-averaged energy, allowing us to compute the twist-averaged energy using only one calculation [27]. Theconnectivity twist averaging scheme reduces the compu-tational cost by a factor of approximately N s , which isvital for our work here. Removing basis set error:-
Running calculationsin a finite basis set of M orbitals incurs an error referredto as basis set incompleteness error. Previous work hasshown that, for a three-dimensional solid, this error de-cays as 1 /M as M → ∞ [18, 28, 36, 37]. In general, tostudy convergence to the thermodynamic limit, we mustfirst eliminate basis set incompleteness error for each elec-tron number N; then, we must find the way that theseenergies asymptote as N → ∞ . Performing this extrapo-lation in both M and N by brute force is impractical. Weinstead derive a correction for basis set incompletenesserror that we may apply before analyzing convergence tothe thermodynamic limit. Full technical details of thisprocedure are given in the Supplemental Material andsummarized below.Following Shepherd and Gr¨uneis [14, 28], we use afixed ratio between the energy cutoff for the virtual andoccupied space for the correlated calculation: f cut = E cut ,M /E cut ,N . This ratio is set to be consistent acrossa set of calculations that share the same density, but dif-fer in electron number N . Then, a correction is derivedfrom a smaller electron number (here, N = 186) extrap-olated to the complete basis set limit, assuming the inde-pendence of the complete basis set and thermodynamiclimits. Numerical comparison to exact results supportan N − / power law:- Figure 1 shows the finite-sizeerror for two different densities over a range of electronnumbers. To isolate the finite-size error, the differencebetween the complete basis set, finite- N CCD energyand the exact thermodynamic limit energy was taken.For r s = 0 .
01 and 0 .
1, RPA energies [17] first derived byGellmann and Brueckner are known to be exact [16], andCCD can be expected to reproduce these numbers. The r s = 0 .
01 calculations required us to use a minimal basis( f cut = √
2) to be able to afford calculations with largerelectron numbers. The remainder used a larger basis set( f cut = 4), which afforded a better residual basis set in-completeness error after correction. Since we are tryingto confirm which power law predominates in the large N N C o rr e l a t i o n F i n i t e S i ze E rr o r E CC D − E e x a c t , ( m H a/ e l ) r s = 0 . r s = 0 . Figure 1. Finite size errors for r s = 0 .
01 and 0 .
1. The electronnumber ranges for each data set are N = 114 to N = 2042 for r s = 0 .
01 and N = 114 to N = 730 for r s = 0 .
1. The dashedline represent free fit exponents (i.e., α in N α ) of − . − . r s = 0 .
01 and r s = 0 .
1, respectively. .
00 0 .
25 0 .
50 0 . N − . . . . . . . N o r m a li ze d F i n i t e S i ze E rr o r , ( E c o rr − E T D L ) / F , ( a r b . un i t s ) .
00 0 . N − . . . . . r s = 1 . r s = 2 . r s = 5 . r s = 10 . r s = 20 . r s = 50 . r s = 100 . Figure 2. A wide range of r s values have correlation energyconvergence to the thermodynamic limit (TDL) which can befit to E corr = E TDL + F N − / . In the left panel, we plot( E corr − E TDL ) /F . The black line is ( E corr − E TDL ) /F = N − / . Different r s values converge at different rates, butsettle down to an N − / trend at large values of N . In theright panel, the same analysis is shown for N − . limit, we plot the data on a log-log scale. Free powerlaw fits to this data yielded exponents of − . − . N − / asymptote [7, 8].Figure 2 compares correlation energy extrapolationsfor both power laws, N − / and N − , for a wide rangeof r s . These data support using an N − / power law.Furthermore, this figure shows the broad applicability ofthis power law, since it is remarkably consistent betweenthe low- and high-density limits. Another surprise is thatthe onset of the power law does not change substantiallywith changing r s over the range of common metal den-sities. In general, we find that reliable extrapolations(i.e., where the error is < N (cid:39)
54 for 5 . < r s < .
0, althoughthis rises sharply with falling r s (e.g., N >
300 electrons for r s = 1). For completeness, the full thermodynamiclimit data we collected can be found in the SupplementalMaterial, where we further compare the N − and N − / extrapolations. Analysis of the transition structure factorshows a cross-over between N − / and N − :- Wehave now shown a significant quantity of numerical datato support an asymptotic N − / power law in the ther-modynamic limit which should reasonably extend to thetotal energy in other wavefunction methods; however,our findings are at odds with a considerable body ofliterature and analysis which support an N − powerlaw [2, 8–11, 14, 25, 29, 38–49]. If this slow asymp-totic convergence is true exclusively for coupled clustertheory, it represents a significant challenge for applyingthis method to real solids; on the other hand, if it ap-plies for all correlation energy calculations, it is a prob-lem for a wide range of previous calculations that usean N − extrapolation. In order to explain this seem-ing contradiction, we need to analyze a property of thecoupled cluster wavefunction called the transition struc-ture factor, which was recently introduced by Gr¨uneisand coworkers [6, 19, 24]. The transition structure factor S G is defined by re-writing the correlation energy ex-pression as: E corr = (cid:80) G t ijab v ijab = (cid:80) G S G V G , where V G is the electron repulsion integral; both S G and V G areFourier components labelled with a momentum G . In the S G → = 0 limit, S G ∼ G can be shown to give rise to the N − relationship derived in the literature [11, 48, 50]. Inour penultimate section we explain the relationship be-tween S G and the true ground-state wavefunction.In order to analyze how the N − / power law emergesand crosses over to N − , we start with fitting a functionalform for S G . A convenient form of the function is: S G ∝ G + λ ) G . (4)Our inspiration for this fit is that it comes from a screenedpotential, ˜ v = e − λr [51], and includes a G term from thesums over the occupied manifold at small G [50]. Thisfunctional form fits the numerically-determined structurefactor for the UEG. If the smallest G present in the cal-culation is given by G (cid:48) , then F ( G (cid:48) ) = (cid:82) G (cid:48) S G V G G dG (where the factor G comes from the G-space volume el-ement in 3D) gives the finite-size error in the correlationenergy. Since G (cid:48) = πL ∝ N − / , the G (cid:48) derivative of F ( G (cid:48) ) gives the slope of the convergence into the TDL ifa N − / power law is assumed. The remaining FSE afterextrapolation is: F ( G (cid:48) ) − G (cid:48) ddG (cid:48) F ( G (cid:48) ) . (5)An analogous expression can be evaluated for an extrap-olation based on G ∝ N − .Figure 3 shows the results of this analysis, taking λ = 1(and ignoring energy prefactors for simplicity, yieldingarbitary units). Two regimes can be seen on either sideof the minimum in S ( G (cid:48) ). For G (cid:48) values to the right ofthe minimum (i.e., when the system size is small), the N − / extrapolation falls with increasing N , and reachesthe thermodynamic limit more quickly than the N − ex-trapolation. However, as G (cid:48) decreases and the systemsize gets larger, the N − / extrapolation eventually over-shoots the true thermodynamic limit. For G (cid:48) to the leftof the minimum (i.e., when the system size is large), the N − extrapolation becomes closer to the thermodynamiclimit, and the N − / extrapolation rises with increasing N . Although N − is clearly the correct extrapolation, atlarge G (cid:48) , it gives larger errors than the N − / extrapola-tion. The way in which the N − / extrapolation behaveswith increasing N tells us which regime we are in. Ifthe prediction from the N − / extrapolation decreaseswith increasing N , we are a long way off the true N − asymptote, and N − / extrapolation gives the best re-sult; however, if the N − / extrapolation increases withincreasing N , then the N − extrapolation gives the bestresult. Practical consequences for calculations:-
Theseanalytical results can also be seen in our real extrap-olation data, shown in Fig. 4. In both r s = 0 . r s = 5 .
0, the N − / extrapolated result first falls, andthen rises again.This shows remarkable visual agreementwith Fig. 3.To evaluate the transferability of this model to realcalculations, we found the G (and thus the N ) corre-sponding to the minimum in the transition structure fac-tor (and thus the cross-over point from the N − / regimeto the N − regime) for each density. For metallic densi-ties between r s = 5 . r s = 20 .
0, the minimum wasfound at N = 100, and the onset of the best N − behav-iors required calculations using 150 < N < s -block metal and thus a close analog tothe UEG, we would expect to need k -point meshes finerthan 5 × × N − extrapolation. Coarser meshes,such as 2 × × × ×
3, are better extrapolatedwith an N − / power law. Relationship to the true ground-statewavefunction:-
The coupled cluster amplitudes, t ijab , are related to the true wavefunction: Whencoupled cluster theory is exact, these are the coeffi-cients in the intermediate-normalized full configurationinteraction (FCI) wavefunction, i.e. t ijab → c ijab where c ijab = (cid:104) Ψ | ˆ a † a ˆ a † b ˆ a i ˆ a j | Ψ HF (cid:105) ; ˆ a are creation andannihilation operators and Ψ HF is the Hartree–Fock(HF) wavefunction. The relationship t ijab = c ijab isvalid in the high density limit; we believe it is goodenough to discuss thermodynamic limit and basis setextrapolations for all densities. Thus, the transitionstructure factor is a property (albeit non-observable) ofthe exact ground-state wavefunction that is a subset of . . . . . . G or G ∝ N − − . . . . . . . F i n i t e S i ze E rr o r ( a r b . un i t s ) N − / N − E corr − E TDL S G − . − . − . − . − . . S G (a.u) Figure 3. An analytical model for the transition structurefactor is used to derive finite-size effects remaining after ex-trapolation to the thermodynamic limit (TDL), in order tocompare the effectiveness of N − / and N − extrapolations.The minimum in the transition structure factor (shown withthe green line), serves as a cross-over point that indicateswhich power law gives a better TDL extrapolation. Symbolicmanipulations were performed in Mathematica [52] .
000 0 .
005 0 . N − − − − − − − C o rr e l a t i o n E n e r g y ( m H a/ e l ) r s = 0 . − / N − .
000 0 .
005 0 . N − − − − − r s = 5 . − / N − . . . . . . Figure 4. Data from r s = 0 . r s = 5 . N , data in the extrap-olated curve (labelled with the power-law) is derived fromconsidering the previous 5 points. Colors match Fig. 3. the transition density matrix between the HF wavefunc-tion and the FCI wavefunction. As such, this propertygoverns convergence behaviors in every high-accuracypost-mean-field method.In a further similarity between FCI and coupled clustertheory, the physical limit of S G → = 0 means that thereis no v G → contribution to the energy, i.e., that the HFMadelung constant does not affect the correlation energy.This suggests a different route toward calculating a cor-rection that would help the convergence of the correlationenergy; we believe that a Madelung-like approach couldbe successful in remedying the error in the correlationenergy where an S G → term is included in truncations of S G . Concluding remarks:-
In summary, we have per-formed coupled cluster calculations on the high-densityuniform electron gas with the intent to derive the powerlaw governing the convergence of energies to the TDL.This was motivated by a recent dilemma that has arisenin the literature around whether an N − / or an N − power law matches convergence of the correlation en-ergy to the thermodynamic limit. By developing andapplying a correction for basis set incompleteness er-ror, complete-basis-set-corrected twist-averaged correla-tion energies show a clear N − / trend in their conver-gence to the thermodynamic limit. These results and the-oretical observations support recent literature that hasshown that this power law is obeyed for insulators andsemiconductors treated in (finite) periodic Gaussian ba-sis sets [7] and the spin-polarized gas [8]. However, whenwe examined the transition structure factor for the UEG,we found that our calculations were distant from the N − asymptote, even at the large system sizes we were ableto reach. We showed analytically that the observation of N − / is a result of using too small an electron numberto resolve the minimum in the transition structure fac-tor (in other words, insufficiently small G to show theappropriate screening in the interaction), and that thepower law reverts to a N − behavior at large N . Theobservation of an N − / power law behavior is relatedto a property derived from the ground-state wavefunc-tion in intermediate normalization and means that thisresult should generalize to correlation energy calculationsacross a broad range of methods.We believe this conclusively explains recent observa-tions related to N − / power laws [7], placing them inthe context of the general expectation that correlation en-ergies, as a component of the total energy, would followan N − power law [2, 8–11, 14, 25, 29, 38–49], (although N − / can also be proposed based on the exchange en-ergy [9, 41]). There are, of course, more sophisticatedmethods to make corrections or models for how the to-tal energy approaches the thermodynamic limit beyondextrapolation [9–11, 25, 38–40, 42, 44, 45, 48]. We arguethat our work is important even in the context of theseapproaches, as it shows that coupled cluster theory hasthe same power law dependence at the doubles level asthese other wavefunction methods.We also discussed that the N − / extrapolation provesto be a powerful predictive tool for whether an N − ex-trapolation would be valid in the absence of analyzingthe transition structure factor directly. This is signifi-cant, because there are a large number of people nowtrying to use coupled cluster theory for materials de-sign [1, 2, 4, 6, 7, 20–22, 24, 53–59]. Finally, we believethis is the first attempt to derive such a simple power lawfor the energy prior to the N − asymptote that also sig-nals which extrapolation regime is being simulated, andthat this general observation is useful well beyond justthe coupled cluster correlation energy. Acknowledgements:-
We gratefully acknowledge theUniversity of Iowa for funding and computer resourcesthrough the University of Iowa Informatics Initiative.We are also grateful to Neil Drummond, MatthewFoulkes, Andreas Gr¨uneis, and Lubos Mitas for theiremails and comments on this work. We thank AliAlavi, David Ceperley, Richard Needs, and Pablo Lopez Rios for useful conversations on related top-ics. Code used throughout this was a locally modifiedversion of a github repository used in previous workhttp://github.com/jamesjshepherd/uegccd [18, 23]. [1] A. Gr¨uneis, Phys. Rev. Lett. , 066402 (2015).[2] G. H. Booth, A. Grneis, G. Kresse, and A. Alavi, Nature , 365 (2013).[3] K. Liao, X.-Z. Li, A. Alavi, and A. Grneis, npj Compu-tational Materials , 1 (2019).[4] A. M. Lewis and T. C. Berkelbach, Physical Review Let-ters , 226402 (2019).[5] X. Wang and T. C. Berkelbach, Journal of Chemical The-ory and Computation , 3095 (2020), publisher: Amer-ican Chemical Society.[6] T. Gruber, K. Liao, T. Tsatsoulis, F. Hummel, andA. Gr¨uneis, Phys. Rev. X , 021043 (2018).[7] J. McClain, Q. Sun, G. K.-L. Chan, and T. C. Berkel-bach, J. Chem. Theory Comput. , 1209 (2017).[8] M. Ruggeri, P. L. R´ıos, and A. Alavi, Phys.Rev. B (2018), 10.1103/PhysRevB.98.161105, arXiv:1804.04938.[9] N. D. Drummond, R. J. Needs, A. Sorouri, and W. M. C.Foulkes, Phys. Rev. B , 125106 (2008).[10] S. Chiesa, D. M. Ceperley, R. M. Martin, and M. Holz-mann, Phys. Rev. Lett. , 076404 (2006).[11] M. Holzmann, R. C. Clay, M. A. Morales, N. M. Tub-man, D. M. Ceperley, and C. Pierleoni, Phys. Rev. B , 035126 (2016).[12] R. F. Bishop and K. H. L¨uhrmann, Phys. Rev. B ,3757 (1978).[13] G. E. Scuseria, T. M. Henderson, and D. C. Sorensen,J. Chem. Phys. , 231101 (2008).[14] J. J. Shepherd and A. Gr¨uneis, Phys. Rev. Lett. ,226401 (2013).[15] J. J. Shepherd, T. M. Henderson, and G. E. Scuse-ria, Range Separated Brueckner Coupled Cluster DoublesTheory , arXiv e-print 1310.6425 (2013).[16] M. Gell-Mann and K. A. Brueckner, Phys. Rev. , 364(1957).[17] L. Onsager, L. Mittag, and M. J. Stephen, Annalen derPhysik , 71 (1966).[18] J. J. Shepherd, T. M. Henderson, and G. E. Scuseria,Phys. Rev. Lett. , 133002 (2014).[19] K. Liao and A. Gr¨uneis, J. Chem. Phys. , 141102(2016).[20] M. Motta, S. Zhang, and G. K.-L. Chan, Physical Re-view B , 045127 (2019), arXiv: 1905.00511.[21] A. Pulkin and G. K.-L. Chan, arXiv:1909.10886 [cond-mat] (2019), arXiv: 1909.10886.[22] Q. Sun, T. C. Berkelbach, J. D. McClain, and G. K.-L.Chan, J. Chem. Phys. , 164119 (2017).[23] J. J. Shepherd, T. M. Henderson, and G. E. Scuseria, J.Chem. Phys. , 124102 (2014).[24] A. Irmler, A. Gallo, F. Hummel, and A. Grneis, PhysicalReview Letters , 156401 (2019).[25] L. M. Fraser, W. M. C. Foulkes, G. Rajagopal, R. J.Needs, S. D. Kenny, and A. J. Williamson, Phys. Rev.B , 1814 (1996).[26] This term represents the exchange energy of an electron interacting with its own periodic image with a neutraliz-ing background. This value is provided for other practi-tioners to reproduce our numbers; v M at all densities andbox lengths can be calculated by re-scaling this value.[27] T. N. Mihm, A. R. McIsaac, and J. J. Shepherd, TheJournal of Chemical Physics , 191101 (2019).[28] J. J. Shepherd, J. Chem. Phys. , 031104 (2016).[29] C. Lin, F. H. Zong, and D. M. Ceperley, Phys. Rev. E , 016702 (2001).[30] T. Gruber and A. Gr¨uneis, Phys. Rev. B , 134108(2018).[31] L. Maschio, D. Usvyat, F. R. Manby, S. Casassa,C. Pisani, and M. Schtz, Phys. Rev. B , 075101 (2007).[32] F. H. Zong, C. Lin, and D. M. Ceperley, Physical ReviewE , 036703 (2002).[33] C. Pierleoni, D. M. Ceperley, and M. Holzmann, PhysicalReview Letters , 146402 (2004).[34] E. Mostaani, N. Drummond, and V. Falko, Physical Re-view Letters , 115501 (2015).[35] S. Azadi and W. M. C. Foulkes, arXiv:1910.06814 [cond-mat, physics:physics] (2019), arXiv: 1910.06814.[36] J. J. Shepherd, G. H. Booth, and A. Alavi, J. Chem.Phys. , 244101 (2012).[37] J. J. Shepherd, A. Gr¨uneis, G. H. Booth, G. Kresse, andA. Alavi, Phys. Rev. B , 035111 (2012).[38] A. J. Williamson, G. Rajagopal, R. J. Needs, L. M.Fraser, W. M. C. Foulkes, Y. Wang, and M.-Y. Chou,Physical Review B , R4851 (1997).[39] P. R. C. Kent, R. Q. Hood, A. J. Williamson, R. J. Needs,W. M. C. Foulkes, and G. Rajagopal, Phys. Rev. B ,1917 (1999).[40] R. Gaudoin and J. M. Pitarke, Physical Review B ,155105 (2007).[41] H. Kwee, S. Zhang, and H. Krakauer, Phys. Rev. Lett. , 126404 (2008).[42] T. Dornheim, S. Groth, T. Sjostrom, F. D. Malone,W. Foulkes, and M. Bonitz, Physical Review Letters , 156403 (2016).[43] D. M. Ceperley and B. J. Alder, Phys. Rev. B , 2092(1987).[44] Y. Kwon, D. M. Ceperley, and R. M. Martin, Phys. Rev.B , 6800 (1998).[45] W. M. C. Foulkes, L. Mitas, R. J. Needs, and G. Ra-jagopal, Rev. Mod. Phys. , 33 (2001).[46] I. G. Gurtubay, R. Gaudoin, and J. M. Pitarke, J. Phys.:Condens. Matter , 065501 (2010).[47] V. S. Filinov, V. E. Fortov, M. Bonitz, and Z. Mold-abekov, Phys. Rev. E , 033108 (2015).[48] M. Holzmann, B. Bernu, and D. M. Ceperley, Journalof Physics: Conference Series , 012020 (2011), pub-lisher: IOP Publishing.[49] D. Ceperley, Phys. Rev. B , 3126 (1978).[50] R. D. Mattuck, A guide to Feynman diagrams in themany-body problem , 2nd ed., Dover books on physics andchemistry (Dover Publications, New York, 1992).[51] A. Gr¨uneis, J. J. Shepherd, A. Alavi, D. P. Tew, andG. H. Booth, J. Chem. Phys. , 084112 (2013).[52] I. Wolfram Research,