Real-frequency response functions at finite temperature
RReal-frequency response functions at finite temperature
I. S. Tupitsyn, A. M. Tsvelik, R. M. Konik, and N. V. Prokof’ev Department of Physics, University of Massachusetts, Amherst, MA 01003, USA Condensed Matter Physics and Materials Science Division,Brookhaven National Laboratory, Upton, NY 11973-5000, USA
Building on a number of previous developments we show that the Diagrammatic Monte Carlotechnique allows one to compute finite temperature response functions directly on the real-frequencyaxis within any field-theoretical formulation of the interacting fermion problem. There are nolimitations on the type and nature of the system’s action or whether partial summation and self-consistent treatment of certain diagram classes are used. By eliminating the need for numericalanalytic continuation from a Matsubara representation, our scheme opens the door for studies ofspectral densities of arbitrary complexity with controlled accuracy. For illustrative purposes weconsider the problem of the plasmon line-width in a homogeneous electron gas (jellium).
Introduction.
The promise of the Diagrammatic MonteCarlo (DiagMC) technique–stochastic sampling of high-order connected Feynman diagrams with extrapolationto the infinite diagram order limit–in solving the com-putational complexity problem for interacting fermions[1] ultimately rests on our ability to formulate a field-theoretical approach with a quickly converging series ex-pansion. As is often the case in the strongly corre-lated regime, an expansion based on the original inter-action potentials, V , and “bare” fermion propagators, G , does not converge. To proceed, the problem is trans-formed identically by incorporating certain classes of di-agrams and interaction effects into an alternative “start-ing point”. This introduces new effective propagators, e G , interactions, U , and counter terms, Λ, in terms ofwhich an alternative diagrammatic expansion is formu-lated. Shifted and homotopic action tools [2, 3] allowone to achieve this goal generically by expressing finalanswers as Taylor series in powers of the auxiliary param-eter ξ , with ξ = 1 corresponding to the original problem,and ensuring that the resulting series converge for any ξ < G and V , the summation over all inter-nal fermionic Matsubara frequencies, ω n = 2 πT ( n + 1 / n , can be performed analytically with thehelp of the Cauchy formula T X n M Y j =1 iω n − a j = M X j =1 n j M Y s = j a j − a s , (1)and automated for arbitrary high-order diagrams. Here n j ≡ n ( a j ) = [ e a j /T +1] − is the Fermi distribution func-tion. [For imaginary time implementation see Ref. [6].]The Wick rotation of external frequency from the imag-inary to the real axis is then performed analytically byreplacing i Ω s with Ω + iη , where η → η are used; ifnecessary, the result is extrapolated to the η → G = ( iω n − (cid:15) k + µ ) − , where (cid:15) k is the bare dispersion relation and µ isthe chemical potential (we suppress the spin index forbrevity), and (ii) a frequency independent interaction po-tential V . These requirements are not satisfied whenthe diagrammatic expansion is performed in terms ofdressed/renormalized propagators and retarded effectiveinteractions to produce convergent series in the stronglycorrelated regimes. Even if e G and U have transparentanalytical structure in the Matsubara representation, thesummation over all internal Matsubara frequencies can-not be performed analytically any more. For example, inthe random phase approximation (RPA) for the homo-geneous electron gas (jellium), the polarization operatorin the effective screened interaction, U − = V − − Π,is approximated by the finite-temperature version of theLindhard function [7]. No diagram with U -lines, as is,can be summed over Matsubara frequencies analytically.It appears that conditions for performing real-frequency simulations are incompatible with the generictools needed to obtain convergent series expansions. Inthis work we present a simple solution to this dilemma a r X i v : . [ c ond - m a t . s t r- e l ] F e b and formulate an approach that allows one to com-pute real-frequency response functions within an arbi-trary field-theoretical setup. To demonstrate how ourapproach works in practice, we compute the plasmonline-width, γ pl , in the jellium model as a function of mo-mentum and temperature. The problem of plasmon de-cay is under active study because of its importance foroptoelectronics, photovoltaics, photocatalysis, and otherapplications (see, for instance, Refs. [8–14] and literaturetherein). Contrary to solid state materials where inter-band transitions and Umklapp processes are possible(and thus the plasmon line-width can be obtained withinthe lowest skeleton order diagrams, in the so-called GWapproximation - see, for instance, Refs. [15, 16]), we findthat meaningful results for γ pl in jellium crucially dependon vertex corrections. Real-frequency finite-temperature technique.
To per-form the Wick rotation by the substitution i Ω s → Ω+ iη ,the function in question has to be known analytically.The key observation leading to solution is that at anypoint in the DiagMC simulation, the propagators and in-teractions used to express the diagram’s contribution areassumed to be known , either analytically or numerically(from relatively simple auxiliary simulations). The firststep is to convert this knowledge into spectral densitiesand use them to express ˜ G and U via˜ G ( k , iω n ) = 1 π Z ∞−∞ du A ( k , u ) iω n − u ; (2) U ( k , iω s ) = V ( k ) + 1 π Z ∞−∞ dv D ( k , v ) iω s − v , (3)with bosonic Matsubara frequencies ω s = 2 πT s . Thesecond step is to rewrite all diagrammatic contributionsin terms of the A and D functions. This will add integra-tions over a set of u and v variables on top of momentum(spin) integrations (sums), which is not a problem forMonte Carlo methods. However, the dependence of theintegrand on Matsubara frequencies is again a product ofsimple poles , meaning that exact summation over all in-ternal Matsubara indexes can be performed analyticallyand the result rotated to the real-frequency axis.The rest of this work is devoted to the explicit demon-stration of how the proposed scheme works in practiceby considering the problem of the plasmon lifetime injellium. Starting point.
First, we need to construct ˜ G and U .The jellium model is defined as the homogenous electrongas on a positive neutralizing background H = X i k i m + X i 0, followingthe standard RPA result in this limit.To complete the setup, we need to compute the spec-tral density D in Eq. (3). For this it is sufficient to knowthe real and imaginary parts of ˜Π = Π (0) + Π (1) = R + iI on the real frequency axis. If we split the total spectraldensity into the electron-hole continuum, D e − h , and thesingular plasmon pole contribution, D pl , then (see alsoRef. [18]) D e − h = − I [( V − − R ) + I ] ; (5) D pl = πr pl ( Q ) δ (Ω − ω pl ( Q )) if I ( Q ) = 0 , (6)where r pl = 1 / | ∂R/∂ Ω | ω pl ( Q ) is the pole residue. Af-ter summation over Matsubara frequencies, the real-frequency result for ˜Π reads:˜Π = − X p F p + Q , p +2 X p , k Y ( p − k ) F p + Q , p F k + Q , k , (7) F q , q = n q − n q Ω − (cid:15) q + (cid:15) q + iη . (8)We evaluated momentum integrals in Eq. (7) by stan-dard Monte Carlo methods on a dense mesh of Q and Ωpoints for several values of κ . The optimized perturba-tion theory strategy [19, 20] would be to choose κ in sucha way that the answer computed up to a given order ofexpansion is least sensitive to its arbitrary value. Pre-vious work in this vein [17] considered static propertiesonly. For the fully dynamic calculation, one is furtherrestricted by the condition that the spectral functionsneed to be positive for any Ω > 0. With respect to thelow-frequency behavior, optimal values of κ would corre-spond to the extrema of the ˜Π(0 , , κ ) curves, shown inlower insets of Fig. 2 for r s = 2 and r s = 4. The factthat both maxima are broad can be used to choose largervalues of κ without loss of accuracy in order to guaranteethat Im ˜Π(Ω > < 0. Indeed, unless κ is large enough,Im ˜Π becomes positive in a finite frequency range, see theupper right inset in Fig. 2. Our strategy then is to chooselarge enough κ as close as possible to the extremum of˜Π(0 , , κ ), leading to κ/k F = 1 . κ/k F = 1 . r s = 2 and r s = 4, respectively. The corresponding realand imaginary parts of ˜Π are presented in Fig. 2. We seethat at moderate values of r s the qualitative behavior re-mains similar to that in the RPA. The smoothing of thesingularities in the polarization operator is a temperatureeffect. Plasmon line-width. In our formulation, the lowest-order polarization diagrams contributing to the finiteplasmon lifetime are shown in Fig. 3. To avoid double-counting, one has to subtract Yukawa potentials fromeffective screened interactions, because the correspond-ing contributions are already included in the definitionsof ˜ G and U functions. The sum of all diagrams in Fig. 3will be denoted as ∆Π.Accounting for the first two diagrams in Fig. 3 wouldbe equivalent to using the so-called GW-approximationperturbatively. It is not surprising then that these twocontributions strongly violate another exact hydrody-namic condition, Π( Q → , Ω = 0) ∝ Q , see Ref. [21], F ReIm (0,0) / RPA (0,0) r S =2 pl F Q k Fpl F ReIm (0,0) / RPA (0,0) r S =4 FIG. 2. (color online) Polarization function, ˜Π = Π (0) + Π (1) ,dependence on frequency at low temperature T /(cid:15) F = 0 . r s = 2, Q/k F = 0 . κ/k F = 1 . 2. Rightpanel: r s = 4, Q/k F = 0 . κ/k F = 1 . 8. Blue and redcurves represent the real and imaginary parts of ˜Π, respec-tively. Lower insets show ˜Π( Q = 0 , Ω = 0) / Π RPA ( Q = 0 , Ω =0) as a function of κ for r s = 2 (left) and r s = 4 (right). In thelimit κ → ∞ Π(0 , 0) saturates at 3 ρ/ ε F equal to 0.05066 inour units. The upper left inset shows the low-momentum partof the plasmon dispersion for r s = 2 within the (i) RPA (solidblack curve), (ii) ˜Π = Π (0) approximation (blue circles), and(iii) ˜Π with vertex correction (red diamonds). (These curveswere extrapolated to η = 0.) The upper right inset shows howIm ˜Π for r s = 4 changes sign for κ/k F < . 8. (These resultsare presented for η/ε F = 0 . (a) (b) (c) U-Y U-Y U-Y FIG. 3. (color online) Lowest-order polarization diagramswithin the formulation based on ˜ G and U . because a similar situation takes place in the GW-approximation [22, 23]. If we were to compute the plas-mon line-width on the basis of the first two diagrams inFig. 3, we would find that the plasmon excitation is com-pletely destroyed at small momenta. Indeed, the datapresented in the left panel of Fig. 4 extrapolate to finitevalues at Q = 0, leading to a divergent contribution aftermultiplication by the Coulomb potential (see also Fig.4of the Supplemental material [24]).It is thus crucial not to miss the vertex correction givenby the diagram (c) in Fig. 3. It compensates diagrams(a) and (b) almost perfectly for all values of Q , andrestores the proper ∝ Q behavior of the (a)+(b)+(c)sum at small momenta, see right panel of Fig. 4. Theinvolved analytical expressions for all diagrams (beforeMonte Carlo integration over internal momenta) can befound in the Supplemental material [24] (see Section I).While their derivation on the basis of Cauchy formula (1)is straightforward, the number of terms rapidly increaseswith the number of frequency dependent lines, not tomention that U functions (3) contain three distinct con-tributions: frequency independent part, plasmon pole,and electron-hole continuum. - I m ( pl ) Q /k F T/ F =0.02 r S = r S = Q /k F Q r S = r S = T/ F =0.02 - I m ( pl ) FIG. 4. (color online) Minus imaginary part of the polariza-tion operator contributions pictured in Fig. 3 as functions ofmomentum for Ω = ω pl ( Q ) and T /ε F = 0 . 02. Left panel:Upper curves are contributions from the sum of diagrams (a)and (b) for r s = 2 (dashed with squares), and r s = 4 (dot-ted with circles). Lower curves are contributions from thediagram (c) for r s = 2 (dashed with triangles), and r s = 4(dotted with diamonds) Right panel: The sum of diagrams(a), (b), and (c) for r s = 2 (red dashed curve with trian-gles) and r s = 4 (blue dashed curve with diamonds). The Q -dependence (black dotted curve) is added for comparison.All error bars are smaller than symbol sizes. After evaluating the imaginary part of ∆Π( Q, Ω = ω pl ( Q )) we obtain the plasmon line-width from γ pl ( Q, T ) = − r pl ( Q, T )Im∆Π(Q , ω pl (Q , T) , T) . (9)[At small momenta r pl ≈ V ω pl / ∝ Q − .] Since the finalresult for γ pl is much smaller than ω pl there is no needfor performing a frequency scan. We have verified thatthe answer does not change when Im ∆Π is computed atfrequencies ω pl ± γ pl .Our final results for the plasmon line-width on thebasis of diagrams with one U -line are discussed inFig. 5. All data are presented as dimensionless ratios γ pl ( Q, T ) /ω pl ( Q, T ) to immediately see when plasmonexcitations remain well-defined. This appears to be thecase all the way to the plasmon spectrum end point forboth values of r s when the temperature is low. The line-width saturates to a finite value in the Q → Q -dependence of Im∆Π is compensated bythe divergence of the Coulomb potential present in thedefinition of the plasmon residue.The answer is also finite in the T → e − h ) excitations that overlaps with the plasmonpeak [18]. Thus there exist kinematically allowed decaychannels for Q = 0 plasmons excited from the groundstate of the system. Somewhat surprising is the fact thatthe line-width remains rather small even for large valesof r s . Finite-temperature corrections are linear at values T /ε F (cid:28) T /ε F > Q /k F r S = F =0.200.10 pl pl pl pl T/ F Q /k F =0.40.30.20.1 T Q /k F pl pl T/ F =0.100.05r S = pl pl T/ F Q /k F =0.60.50.40.30.1 T FIG. 5. (color online) Plasmon line-width to plasmon fre-quency ratio as a function of momentum at different temper-atures for r s = 2 (left panel) and 4 (right panel). The blackdotted lines in both panels show results extrapolated to the T → Q is shown in insets. Higher order contributions. Our reformulation of thediagrammatic expansion in terms of ˜ G and U is exact,and one can proceed with computing higher-order dia-grams using standard rules. We illustrate some of thesecond-order diagrams and process them in the Supple-mental material [24] (see Section III). While summa-tion over internal Matsubara frequencies allows one toperform calculations directly on the real-frequency axis,it also brings additional computational challenges. Us-ing two next-order diagrams as an example, we demon-strate in the Supplemental material that processing Mat-subara sums “by hand” quickly leads to expressions ofoverwhelming complexity (the remaining momentum in-tegrals are done by standard Monte Carlo techniques).Since the Cauchy formula (1) is recursive, it should bepossible to fully automate the process, similarly to whatwas done in Ref. [5] for the case when only fermionicpropagators were frequency dependent.Given that certain groups of diagrams feature strongcompensation, an efficient algorithm would need to com-bine them analytically (see also Ref. [17]). We see clearadvantages in implementing the recursive scheme for ob-taining Taylor expansions from skeleton diagrams [25],because it automatically groups irreducible diagrams andreduces the number propagators and interaction lines.This scheme also significantly simplifies processing ofcounter terms, and eliminates higher order poles in Mut-subara sums. Conclusion. We report a simple solution to the prob-lem of computing finite-temperature response functionson the real frequency axis using Feynman diagrams thateliminates the need for infamous numerical analyticalcontinuation. Our solution works for arbitrary field-theoretical formulation of the problem including dressed,renormalized, and self-consistent treatments required forproducing convergent expansions. Spectral densities (ofarbitrary complexity) for experimentally relevant observ-ables can be computed with accuracy that was neverpossible before, opening the door for numerous appli-cations involving optical conductivity, resonant inelasticX-ray spectroscopy, neutron scattering, excitation life-times, etc.To illustrate how the technique works, we used it tocompute the leading processes contributing to the finiteplasmon line-with within the jellium model, and stud-ied the line-width dependence on momentum and tem-perature for moderate values of the Coulomb r s param-eter. The increase of the interaction strength leads toa decrease of the plasmon life time, but nevertheless theplasmon remains well defined. One important qualitativeresult is the necessity to include vertex corrections in or-der to ensure the obtained results do not violate generalprinciples.Future work will aim at developing efficient schemes forgenerating and processing real-frequency expressions forhigh-order diagrams to gain full control over systematicerrors resulting from the series truncation. Acknowledgements. AMT, RMK, and IST thank sup-port from the Office of Basic Energy Sciences, Mate-rial Sciences and Engineering Division, U.S. Departmentof Energy under Contract No. DE-SC0012704. NVPthanks support from the Simons Collaboration on theMany Electron Problem. The authors thank JamesLeBlank and Kun Chen for sharing details on their meth-ods and helpful discussions. [1] R. Rossi, N. Prokof’ev, B. Svistunov, K. Van Houcke,and F. Werner, Euro Phys. Lett. , 10004 (2017).[2] R. Rossi, F. Werner, N. Prokof’ev, B. Svistunov, Phys.Rev. B 93 , 161102(R) (2016).[3] A. J. Kim, N. V. Prokof’ev, B. V. Svistunov, and E. Kozik, arXiv:2010.05301[4] O. Goulko, A.S. Mishchenko, L. Pollet, N. Prokof’ev, andB. Svistunov, Phys. Rev. B 95 , 014102 (2017).[5] A. Taheridehkordi, S.H. Curnoe, and J.P.F. LeBlanc,Phys. Rev. B 99 , 035120 (2019); ibid. arXiv preprintarXiv:2004.11091.[6] J. Viˇciˇcevi´c, P. Stipsi´c, and M. Ferrero,arXiv:2011.08226.[7] J. Lindhard, Mat. Fys. Medd. K. Dan. Vidensk. Selsk,28, 1 (1954).[8] J. M. Luther and J. L. Blackburn, Nat. Photonics , 675(2013).[9] C. Clavero, Nat. Photonics , 95 (2014).[10] H. A. Atwater and A. Polman, Nat. Mater. , 205 (2010).[11] S. Linic, P. Christopher, and D. B. Ingram, Nat. Mater. , 911 (2011).[12] S. Mukherjee, F. Libisch, N. Large, O. Neumann, L. V.Brown, J. Cheng, J. B. Lassiter, E. A. Carter, P. Nord-lander, and N. J. Halas, Nano Lett. , 240 (2013).[13] W. Li and J.G. Valentine, Nanophotonics (1), 177,(2017).[14] K. Kolwas, Plasmonics , 1629 (2019).[15] M. Bernardi, J. Mustafa, J.B. Neaton and S.G. Louie,Nat. Commun. , 7044 (2015).[16] R. Sundararaman, P. Narang, A.S. Jermyn, W.A. God-dard III and Harry A. Atwater, Nat. Commun. , 5788(2014).[17] K. Chen and K. Haule, Nat. Commun. , 3725 (2019);ibid. arXiv:2012.03146.[18] I.S. Tupitsyn, A.M. Tsvelik, R.M. Konik, N.V, Prokof’ev,Phys. Rev. B , 102 (7), 075140 (2020).[19] P.M. Stevenson, Phys. Rev. D 23 , 2916 (1981).[20] R.P. Feynman and H.E. Kleinert, Phys. Rev. A 34 , 5080(1986).[21] P. Nozieres and D. Pines, Theory Of Quantum Liquids ,Chapters 2 and 3, Westview Press (1999).[22] B. Holm and U. von Barth, Phys. Rev. B 57 , 2108 (1998).[23] K.V. Houcke, I.S. Tupitsyn, A.S. Mishchenko, N.V.Prokof’ev, Phys. Rev. B 95 (19), 195131 (2017).[24] Supplemental material[25] W. Wu, M. Ferrero, A. Georges, and E. Kozik, Phys.Rev. B 96 , 041105 (2017) upplemental material for ”Real-frequency response functions at finite temperature” I. S. Tupitsyn, A. M. Tsvelik, R. M. Konik, and N. V. Prokof’ev I. ANALYSIS OF LEADING ORDER DIAGRAMS IN FIG. 3 OF MAIN TEXT Let us start with specifying all relevant external and internal variables for the diagram shown in Fig.3(a) of the main text, seeFig. 1. , s Q i , s Q i , n s k Q i i , n m k q i i , n k i , n k i , m q i FIG. 1: (color online) The same polarization diagram as in Fig.3(a) of the main text with all frequency and momentum variables specified. Q and Ω are considered to be external parameters. Spin indices are not shown for brevity. Before we proceed, let us introduce a convenient set of short-hand notations: (cid:15) k , p = (cid:15) ( k ) − (cid:15) ( p ) , n u = n ( u ) = 1 e u/T + 1 , δn k , p = n (cid:15) ( k ) − n (cid:15) ( p ) , n k = dn ( x ) dx | x = (cid:15) ( k ) .N u = N ( u ) = 1 e u/T − , N k , p = N ( (cid:15) k , p ) . Since the interaction line is given by U − Y = V ( q ) − Y ( q ) + (1 /π ) Z ∞−∞ du D ( q , u ) iω m − u , We thus define their contributions to ∆Π ( a ) via ∆Π ( a ) = ∆Π ( a ) V Y + ∆Π ( a ) D . A recursive application of the Cauchy formula forthe V C ( q ) − Y ( q ) term results in ∆Π ( a ) V Y = − X k , q πe κ q ( q + κ ) n k − q − i Ω s − (cid:15) k − Q , k (cid:20) n k − δn k , k − Q − i Ω s − (cid:15) k − Q , k (cid:21) . (1)For the second term we get: ∆Π ( a ) D = − π X k , q Z ∞−∞ duD ( q , u ) ( n k − q − − N u )( i Ω s + (cid:15) k − Q , k )( (cid:15) k − q , k + u ) (2) × (cid:20) n k + δn k , k − Q i Ω s + (cid:15) k − Q , k + n k − n ( (cid:15) k − q + u ) (cid:15) k − q , k + u + n k − Q − n ( (cid:15) k − q + u ) i Ω s − (cid:15) k − q , k − Q − u (cid:21) . (3)Wick rotation to the real-frequency domain is now performed by a simple substitution i Ω s → Ω + iη . For diagram (b) in Fig.3of the main text the analytical expression for ∆Π ( b ) is obtained from that of ∆Π ( a ) by the replacement Ω + iη → − Ω − iη . Fordiagram (c) the final result for ∆Π ( c ) = ∆Π ( c ) V Y + ∆Π ( c ) D reads ∆Π ( c ) V Y = − X k , q πe κ q ( q + κ ) δn k , k − Q ( − Ω − iη − (cid:15) k − Q , k ) δn k − q , k − q − Q ( − Ω − iη − (cid:15) k − q − Q , k − q ) , (4) a r X i v : . [ c ond - m a t . s t r- e l ] F e b ∆Π ( c ) D = − π X k , q Z ∞−∞ du D ( q , u )( − Ω − iη − (cid:15) k − Q , k )( − Ω − iη − (cid:15) k − q − Q , k − q ) (5) × (cid:26) ( n k − q − − N u ) (cid:20) n k − Q − n ( (cid:15) k − q + u ) − Ω − iη + (cid:15) k − q , k − Q + u − n k − n ( (cid:15) k − q + u ) (cid:15) k − q , k + u (cid:21) − ( n k − q − Q − − N u ) (cid:20) n k − Q − n ( (cid:15) k − q − Q + u ) (cid:15) k − q − Q , k − Q + u − n k − n ( (cid:15) k − q − Q + u )Ω + iη + (cid:15) k − q − Q , k + u (cid:21)(cid:27) . II. SECOND ORDER DIAGRAMS BASED ON YUKAWA POTENTIAL AND COUNTER TERMS One may wonder whether having a “built in” plasmon excitation and more elaborate field-theoretical formulation has a sig-nificant effect on the final result. We can answer in the affirmative here. To see this, we consider effective interactions based onthe same Yukawa potential as in the Hartree-Fock treatment of ˜ G (this choice was made in Ref. 1 which provides further detailsabout the expansion) and compute contributions from diagrams having the same structure as diagrams in Fig. 3 of the main text.These diagrams are obtained by replacing U − Y lines with the bubble diagram and the associated counter-term, see Fig. 2. (a) (b) (c) FIG. 2: Upper row: Second-order (with respect to the number of interaction lines) diagrams for the polarization operator with self-energyinsertions on the upper (a) and lower (b) propagator lines and vertex correction (c). The counter-terms are shown in the lower row. Dashedlines represent the Yukawa potential. , s Q i , s Q i , n q i i , k q i , n k i , n k i , r p i , r n p q i i i , n s k Q i i , n q i i , s Q i , s Q i , n s k Q i i , k q i , n k i , n k i , n q i i , n q i i /4 FIG. 3: (color online) Second-order (with respect to the number of interaction lines) diagram for polarization operator with self-energy insertion(left); the counter-term is shown to the right. Dashed lines represent the Yukawa potential. Fig. 3 here provides all details for the diagrams in Fig. 2(a). The counter-term (CT) contribution is the easiest one: Π ( a ) CT = 2 κ πe X k , q Y ( q ) n k − q − i Ω s − (cid:15) k − Q , k (cid:20) n k − δn k , k − Q − i Ω s − (cid:15) k − Q , k (cid:21) . (6)The main diagram has two additional poles and an extra momentum integration: Π ( a ) MD = − X k , p , q Y ( q ) ( n k − q + n p + q − n p + N ( (cid:15) p + q + (cid:15) k − q ))( − i Ω s − (cid:15) k − Q , k )( (cid:15) p + q , p + (cid:15) k − q , k ) (7) × (cid:20) n ( (cid:15) p + q + (cid:15) k − q ) − n k (cid:15) p + q , p + (cid:15) k − q , k − n k + δn k , k − Q − i Ω s − (cid:15) k − Q , k − n ( (cid:15) p + q , p + (cid:15) k − q ) − n k − Q − i Ω s + (cid:15) p + q , p + (cid:15) k − q , k − Q (cid:21) . (8)The complementary expression for diagrams in Fig. 2(b) is obtained from (8) after Wick rotation by the replacement Ω + iη →− Ω − iη . Diagrams in Fig. 2(c) are computed following the same protocol.In Fig. 4 we present results based on calculations of diagrams (a), (b), and (c) in Fig. 2 at r s = 2 and T /ε F = 0 . for Ω = ω pl ( Q ) . Taken separately, these diagrams are physically meaningless because they saturate to a constant value at Q → , and have opposite sign. Only the sum of all three diagrams leads to the proper Q behavior at small Q . Quantitatively,these diagrams predict the plasmon line-width that is nearly a factor of two smaller than results presented in the main text,demonstrating the importance of constructing an expansion with “built in” collective excitations. - I m Q /k F r S =2 , T/ F =0.02 - I m Q /k F Q r S =2, T/ F =0.02 FIG. 4: (color online) Left panel: Minus imaginary part of partial contributions to the polarization operator from diagrams in Fig. 2 as afunctions of momentum at r s = 2 , T /ε F = 0 . , and Ω = ω pl ( Q ) . The sum of diagrams (a) and (b) is shown by red squares. The vertexcorrection (c) is shown by red up-triangles. Right panel: The sum of all three diagrams (red up-triangles) is compared to the Q fit (dottedblack curve). III. SECOND-ORDER DIAGRAMS WITH U POTENTIALS In Fig. 5 we list some of the second-order diagrams that appear in the expansion based on ˜ G and U functions. In whatfollows we demonstrate that higher-order real-frequency calculations for expansions based on retarded effective interactionsare feasible. In this work we deal with all Matsubara sums by hand. The reader can immediately see that the complexity ofexpressions quickly becomes overwhelming. With additional algorithmic developments, this task can be fully automated.As an example, we compute the first two diagrams shown in Fig. 5. Their complete characterization is given in Fig. 6. Firstwe performs sums over indexes n and k to obtain contributions associated with the left and right loops. For the left loops wehave f L ( k ) = 1 i Ω s − (cid:15) k + Q , k (cid:20) δn k , k + q iω m − (cid:15) k + q , k − δn k + Q , k + q iω m − i Ω s − (cid:15) k + q , k + Q (cid:21) . (9)The right loop for diagram (a) makes an identical contribution with k → p . For the right loop in diagram (b) we get f R ( p ) = 1 − i Ω s − (cid:15) p + Q , p (cid:20) δn p , p + q − iω m − (cid:15) p + q , p − δn p + Q , p + q − iω m + i Ω s − (cid:15) p + q , p + Q (cid:21) . (10) (a) (b)(c) (d) (f)(e) FIG. 5: (color online) Examples of second order polarization diagrams in the new formulation. To differentiate between the U − Y functionsthat take care of double-counting and U , we used single-wavy lines for the former in the main text and double-wavy lines for the latter here. , s Q i , s Q i , n k i , n s k Q i i , n m k q i , m s q Q i i , m q i , k s p Q i i , k p i , k m p q i (a) ё , s Q i , s Q i , n k i , n s k Q i i , n m k q i , m s q Q i i , m q i , k s p Q i i , k p i , k m p q i (b) FIG. 6: (color online) Second order polarization diagrams with two fermionic loops. For convenience, in the diagram ( b ) we flip the sign ofthe momentum variable for the right loop and use reflection symmetry (cid:15) ( − p ) = (cid:15) ( p ) = (cid:15) ( p ) . Note that ω m and Ω s are bosonic Matsubarafrequencies. Given that each U function has two dissimilar terms, there are four contributions to the polarization operator: Π (2) a + b = − P k , p , q [ C + C + C + C ] , where C = V C ( Q − q ) V C ( q ) T X m F ; (11) C = 1 π Z ∞−∞ duV C ( Q − q ) D ( q , u ) T X m F iω m − u ; (12) C = 1 π Z ∞−∞ duD ( Q − q , u ) V C ( q ) T X m F− iω m + i Ω s − u ; (13) C = 1 π Z ∞−∞ du Z ∞−∞ dvD ( Q − q , u ) D ( q , v ) T X m F ( − iω m + i Ω s − u )( iω m − v ) ; (14) F = f L ( k , iω m )[ f L ( p , iω m ) + f R ( p , iω m )] . (15)Sums over m lead to functions S , S , S , and S described below. They have to be substituted into Eqs. (11-14) and theresult further integrated over internal momenta to obtain Π (2) a + b . S . For diagram (a), S is given by the product of [( i Ω s − (cid:15) k + Q , k )( i Ω s − (cid:15) p + Q , p )] − and (cid:26) δn k , k + q δn p , p + q [ N p + q , p − N k + q , k ] (cid:15) k + q , k − (cid:15) p + q , p + δn k + Q , k + q δn p + Q , p + q [ N p + q , p + Q − N k + q , k + Q ] (cid:15) k + q , k + Q − (cid:15) p + q , p + Q − δn k , k + q δn p + Q , p + q [ N p + q , p + Q − N k + q , k ] − i Ω s + (cid:15) k + q , k − (cid:15) p + q , p + Q − δn k + Q , k + q δn p , p + q [ N p + q , p − N k + q , k + Q ] i Ω s + (cid:15) k + q , k + Q − (cid:15) p + q , p (cid:27) . (16)Similarly, for diagram (b), S is given by the product of [( i Ω s − (cid:15) k + Q , k )( − i Ω s − (cid:15) p + Q , p )] − and (cid:26) δn k , k + q δn p , p + q [1 + N p + q , p + N k + q , k ] (cid:15) k + q , k + (cid:15) p + q , p + δn k + Q , k + q δn p + Q , p + q [1 + N p + q , p + Q + N k + q , k + Q ] (cid:15) k + q , k + Q + (cid:15) p + q , p + Q − δn k , k + q δn p + Q , p + q [1 + N p + q , p + Q + N k + q , k ] − i Ω s + (cid:15) k + q , k + (cid:15) p + q , p + Q − δn k + Q , k + q δn p , p + q [1 + N p + q , p + N k + q , k + Q ] i Ω s + (cid:15) k + q , k + Q + (cid:15) p + q , p (cid:27) . (17) S . For diagram ( a ) , S is given by the product of [( i Ω s − (cid:15) k + Q , k )( i Ω s − (cid:15) p + Q , p )] − and (cid:26) δn k , k + q δn p , p + q (cid:15) k + q , k − (cid:15) p + q , p × (cid:20) N k + q , k − N u u − (cid:15) k + q , k − N p + q , p − N u u − (cid:15) p + q , p (cid:21) + δn k + Q , k + q δn p + Q , p + q (cid:15) k + q , k + Q − (cid:15) p + q , p + Q × (cid:20) N k + q , k + Q − N u − i Ω s + u − (cid:15) k + q , k + Q − N p + q , p + Q − N u − i Ω s + u − (cid:15) p + q , p + Q (cid:21) − δn k , k + q δn p + Q , p + q − i Ω s + (cid:15) k + q , k − (cid:15) p + q , p + Q × (cid:20) N k + q , k − N u u − (cid:15) k + q , k − N p + q , p + Q − N u − i Ω s + u − (cid:15) p + q , p + Q (cid:21) − δn k + Q , k + q δn p , p + q i Ω s + (cid:15) k + q , k + Q − (cid:15) p + q , p × (cid:20) N k + q , k + Q − N u − i Ω s + u − (cid:15) k + q , k + Q − N p + q , p − N u u − (cid:15) p + q , p (cid:21) (cid:27) . (18)For diagram ( b ) , S is given by the product of [( i Ω s − (cid:15) k + Q , k )( − i Ω s − (cid:15) p + Q , p )] − and (cid:26) δn k , k + q δn p , p + q (cid:15) k + q , k + (cid:15) p + q , p × (cid:20) N u − N k + q , k u − (cid:15) k + q , k − N u + N p + q , p u + (cid:15) p + q , p (cid:21) + δn k + Q , k + q δn p + Q , p + q (cid:15) k + q , k + Q + (cid:15) p + q , p + Q × (cid:20) N u − N k + q , k + Q − i Ω s + u − (cid:15) k + q , k + Q − N u + N p + q , p + Q − i Ω s + u + (cid:15) p + q , p + Q (cid:21) − δn k , k + q δn p + Q , p + q − i Ω s + (cid:15) k + q , k + (cid:15) p + q , p + Q × (cid:20) N u − N k + q , k u − (cid:15) k + q , k − N u + N p + q , p + Q − i Ω s + u + (cid:15) p + q , p + Q (cid:21) − δn k + Q , k + q δn p , p + q i Ω s + (cid:15) k + q , k + Q + (cid:15) p + q , p × (cid:20) N u − N k + q , k + Q − i Ω s + u − (cid:15) k + q , k + Q − N u + N p + q , p u + (cid:15) p + q , p (cid:21) (cid:27) . (19) S . For diagram ( a ) , S is given by the product of [( i Ω s − (cid:15) k + Q , k )( i Ω s − (cid:15) p + Q , p )] − and (cid:26) δn k , k + q δn p , p + q (cid:15) k + q , k − (cid:15) p + q , p × (cid:20) N u + N k + q , k − i Ω s + u + (cid:15) k + q , k − N u + N p + q , p − i Ω s + u + (cid:15) p + q , p (cid:21) + δn k + Q , k + q δn p + Q , p + q (cid:15) k + q , k + Q − (cid:15) p + q , p + Q × (cid:20) N u + N k + q , k + Q u + (cid:15) k + q , k + Q − N u + N p + q , p + Q u + (cid:15) p + q , p + Q (cid:21) − δn k , k + q δn p + Q , p + q − i Ω s + (cid:15) k + q , k − (cid:15) p + q , p + Q × (cid:20) N u + N k + q , k − i Ω s + u + (cid:15) k + q , k − N u + N p + q , p + Q u + (cid:15) p + q , p + Q (cid:21) − δn k + Q , k + q δn p , p + q i Ω s + (cid:15) k + q , k + Q − (cid:15) p + q , p × (cid:20) N u + N k + q , k + Q u + (cid:15) k + q , k + Q − N u + N p + q , p − i Ω s + u + (cid:15) p + q , p (cid:21) (cid:27) . (20)For diagram ( b ) , S is given by the product of [( i Ω s − (cid:15) k + Q , k )( − i Ω s − (cid:15) p + Q , p )] − and (cid:26) δn k , k + q δn p , p + q (cid:15) k + q , k + (cid:15) p + q , p × (cid:20) − N u + N k + q , k − i Ω s + u + (cid:15) k + q , k + N u − N p + q , p − i Ω s + u − (cid:15) p + q , p (cid:21) + δn k + Q , k + q δn p + Q , p + q (cid:15) k + q , k + Q + (cid:15) p + q , p + Q × (cid:20) − N u + N k + q , k + Q u + (cid:15) k + q , k + Q + N u − N p + q , p + Q u − (cid:15) p + q , p + Q (cid:21) − δn k , k + q δn p + Q , p + q − i Ω s + (cid:15) k + q , k + (cid:15) p + q , p + Q × (cid:20) − N u + N k + q , k − i Ω s + u + (cid:15) k + q , k + N u − N p + q , p + Q u − (cid:15) p + q , p + Q (cid:21) − δn k + Q , k + q δn p , p + q i Ω s + (cid:15) k + q , k + Q + (cid:15) p + q , p × (cid:20) − N u + N k + q , k + Q u + (cid:15) k + q , k + Q + N u − N p + q , p − i Ω s + u − (cid:15) p + q , p (cid:21) (cid:27) . (21) S . For diagram ( a ) , S is given by the product of [( i Ω s − (cid:15) k + Q , k )( i Ω s − (cid:15) p + Q , p )( i Ω s − u − v )] − and (cid:26) δn k , k + q δn p , p + q (cid:15) k + q , k − (cid:15) p + q , p × (cid:20) N u + N k + q , k − i Ω s + u + (cid:15) k + q , k − N u + N p + q , p − i Ω s + u + (cid:15) p + q , p + N k + q , k − N v v − (cid:15) k + q , k − N p + q , p − N v v − (cid:15) p + q , p (cid:21) + δn k + Q , k + q δn p + Q , p + q (cid:15) k + q , k + Q − (cid:15) p + q , p + Q × (cid:20) N u + N k + q , k + Q u + (cid:15) k + q , k + Q − N u + N p + q , p + Q u + (cid:15) p + q , p + Q + N k + q , k + Q − N v − i Ω s + v − (cid:15) k + q , k + Q − N p + q , p + Q − N v − i Ω s + v − (cid:15) p + q , p + Q (cid:21) − δn k , k + q δn p + Q , p + q − i Ω s + (cid:15) k + q , k − (cid:15) p + q , p + Q × (cid:20) N u + N k + q , k − i Ω s + u + (cid:15) k + q , k − N u + N p + q , p + Q u + (cid:15) p + q , p + Q + N k + q , k − N v v − (cid:15) k + q , k − N p + q , p + Q − N v − i Ω s + v − (cid:15) p + q , p + Q (cid:21) − δn k + Q , k + q δn p , p + q i Ω s + (cid:15) k + q , k + Q − (cid:15) p + q , p × (cid:20) N u + N k + q , k + Q u + (cid:15) k + q , k + Q − N u + N p + q , p − i Ω s + u + (cid:15) p + q , p + N k + q , k + Q − N v − i Ω s + v − (cid:15) k + q , k + Q − N p + q , p − N v v − (cid:15) p + q , p (cid:21) (cid:27) . (22)For diagram ( b ) , S is given by the product of [( i Ω s − (cid:15) k + Q , k )( − i Ω s − (cid:15) p + Q , p )( i Ω s − u − v )] − and (cid:26) δn k , k + q δn p , p + q (cid:15) k + q , k + (cid:15) p + q , p × (cid:20) − N u + N k + q , k − i Ω s + u + (cid:15) k + q , k + N u − N p + q , p − i Ω s + u − (cid:15) p + q , p + N v − N k + q , k v − (cid:15) k + q , k − N v + N p + q , p v + (cid:15) p + q , p (cid:21) + δn k + Q , k + q δn p + Q , p + q (cid:15) k + q , k + Q + (cid:15) p + q , p + Q × (cid:20) − N u + N k + q , k + Q u + (cid:15) k + q , k + Q + N u − N p + q , p + Q u − (cid:15) p + q , p + Q + N v − N k + q , k + Q − i Ω s + v − (cid:15) k + q , k + Q − N v + N p + q , p + Q − i Ω s + v + (cid:15) p + q , p + Q (cid:21) − δn k , k + q δn p + Q , p + q − i Ω s + (cid:15) k + q , k + (cid:15) p + q , p + Q × (cid:20) − N u + N k + q , k − i Ω s + u + (cid:15) k + q , k + N u − N p + q , p + Q u − (cid:15) p + q , p + Q + N v − N k + q , k v − (cid:15) k + q , k − N v + N p + q , p + Q − i Ω s + v + (cid:15) p + q , p + Q (cid:21) − δn k + Q , k + q δn p , p + q i Ω s + (cid:15) k + q , k + Q + (cid:15) p + q , p × (cid:20) − N u + N k + q , k + Q u + (cid:15) k + q , k + Q + N u − N p + q , p − i Ω s + u − (cid:15) p + q , p + N v − N k + q , k + Q − i Ω s + v − (cid:15) k + q , k + Q − N v + N p + q , p v + (cid:15) p + q , p (cid:21) (cid:27) . (23)To move to real frequencies, replace i Ω s with Ω + iη before all integrations. In Fig. 7 we show our results for − Im Π (2) a + b ( Q , Ω = ω pl ( Q )) for r s = 2 for T /ε F = 0 . and . . All by itself, this partial contribution is not physicalbecause it has a wrong sign (the same was true for lower-order vertex corrections). K. Chen and K. Haule, Nat. Commun. , 3725 (2019); ; ibid. arXiv:2012.03146. Q /k F T/ F =0.02T/ F =0.2 - Im a+b r S =2 FIG. 7: (color online) − Im Π a + b as a function of momentum Q for r s = 2 and T /ε F = 0 . and .2