On the Chemical Potential of Many-Body Perturbation Theory in Extended Systems

Finite-temperature many-body perturbation theory in the grand-canonical ensemble is fundamental to numerous methods for computing electronic properties at nonzero temperature, such as finite-temperature coupled-cluster. In most applications it is the average number of electrons that is known rather than the chemical potential. Expensive correlation calculations must be repeated iteratively in search for the interacting chemical potential that yields the desired average number of electrons. In extended systems with mobile charges the situation is particular, however. Long-ranged electrostatic forces drive the charges such that the average ratio of negative and positive charges is one for any finite chemical potential. All properties per electron are expected to be virtually independent of the chemical potential, as they are in an electric wire at different voltage potentials. This work shows that per electron, the exchange-correlation free energy and the exchange-correlation grand potential indeed agree in the infinite-size limit. Thus, only one expensive correlation calculation suffices for each system size, sparing the search for the interacting chemical potential. This work also demonstrates the importance of regularizing the Coulomb interaction such that each electron on average interacts only with as many electrons as there are electrons in the simulation, avoiding interactions with periodic images. Numerical calculations of the warm uniform electron gas have been conducted with the Spencer–Alavi regularization employing the finite-temperature Hartree approximation for the self-consistent field and linearized finite-temperature direct-ring coupled-cluster doubles for treating correlation.


BACKGROUND
In the warm-dense matter (WDM) regime a large number of many-body states are thermalized. Also, the density is sufficiently large to require a quantum mechanical treatment of the electrons interacting with each other. WDM conditions are found, for instance, during inertial confinement fusion (ICF), in the metallic phase of hydrogen in gas giants, or in matter interacting with high intensity laser fields. 1 Even at room temperature the thermal energy must be considered to be large compared to the vanishing band gap of bulk metals.
The mobility of electrons at warm-dense conditions poses challenges for ab initio simulations of extended systems that are absent in zero-temperature calculations. Unlike at zero temperature, the number of electrons in a volume of fixed shape fluctuates. Periodic boundary conditions, usually employed in extended systems, cannot model such fluctuations at length scales beyond the size of the simulation cell. In reality, the charges would move from one cell to the neighboring cell, keeping the average charge constant. Under periodic boundary conditions, on the other hand, charges can only appear or disappear simultaneously in all periodic images of the simulation cell. However, configurations of periodically repeating netcharged cells are not allowed due to the diverging electrostatic energy per volume. There are mainly two methods in current state-of-the-art ab initio simulations at warm-dense conditions to circumvent this difficulty: (i) The simulation is done in the canonical ensemble where electrons are not permitted to enter or leave the simulated volume. While this ensures charge neutrality it also reduces the number of possible configurations, affecting the system's entropy. 2 Path-integral quantum Monte Carlo (PIQMC) calculations are usually conducted in the canonical ensemble. 3,4 (ii) Another possibility is to disregard the parts of the electrostatic interaction stemming from the average electron and background densities, thus removing the divergence. This allows for grand-canonical simulations with a fluctuating number of electrons including its effect on the entropy. Perturbation-theory-related calculations usually apply this method 5,6 following the work of Kohn and Luttinger, in particular the assumption for arriving at eq 20 in ref 7. A physical justification for this procedure would be if the fluctuations of the positive background were fully correlated with the fluctuations of the electrons. Different mobilities of electrons and ions, however, question this assumption.
In this work, a third alternative is studied to treat long-range electrostatic interactions under periodic boundary conditions. Liang and co-workers 8 have studied classical simulations of mobile electrostatically interacting particles under periodic boundary conditions. They look at the pair correlation function and observe the theoretically expected Debey−Hueckel screening at long distances only under two conditions: (i) when simulating in the grand-canonical ensemble, and (ii) when limiting the range of the electrostatic interaction, such that the particles do not interact with all of their own periodic images. Here, the Coulomb interaction is truncated spherically such that the sphere's volume agrees with the volume of the simulation cell. This modification is called regularization of the Coulomb interaction and it turns divergent terms into finite terms. Care must be taken that the infinite-size limits of the computed quantities are not dependent on the details of the regularization scheme, as discussed in section 2.
A spherical truncation scheme has already been developed by Spencer and Alavi 9 to avoid spurious Fock-exchange interactions of the electrons with their periodic images for zerotemperature calculations as an alternative to other methods treating the occurring integrable singularity. 10,11 This work applies the truncation scheme to all parts of the electrostatic interaction in the self-consistent field calculations, as well as in the subsequent correlation calculation. Other regularization schemes that limit the interaction range are also possible, such as the minimal image convention (MIC) for atom centered orbitals, or the Wigner−Seitz truncation scheme. 12, 13 For pointlike charges the spherical truncation is not continuous which may pose difficulties when considering different atomic configurations.
Related Work. Finite-temperature many-body perturbation theory (FT-MBPT) offers an elementary framework for ab initio calculations of WDM. 5,6,14−16 Numerous approximation schemes employ thermal MBPT, such as thermal secondorder MBPT, 17−20 finite-temperature random phase approximation, 21−24 Green's function based methods, 25,26 as well as some finite-temperature generalizations of coupled-cluster methods. 27−30 An alternative formulation of the coupled-cluster methods has been brought forward recently in the framework of thermo-field dynamics. 31−33 Finite-temperature perturbation theory is originally formulated in the grand-canonical ensemble; however, formulations in the canonical ensemble exist. 34, 35 Equally, thermo-field dynamics can be employed in the canonical ensemble. 36 Analogous to ab initio calculations at zero-temperature, thermal Hartree−Fock and density functional theory (DFT) calculations are among the most widely used methods serving also as self-consistent field (SCF) reference for FT-MBPT. 37−39 For DFT it is in general not sufficient to use a zero-temperature exchange-correlation functional and introduce temperature merely by smearing. Temperature must be a parameter of the exchange-correlation functional. 40 At higher temperatures, a large number of one-body states is occupied with non-negligible probabilities. Orbital-free density functional theories (ofDFT) aim at mitigating this with functionals that do not depend on the usual Kohn−Sham orbital description of DFT. 41,42 Canonical or grand-canonical full configuration interaction methods can be used for benchmarking more approximate theories. 43 Finally, various forms of Monte Carlo methods approach the finitetemperature many-body problem in complementary ways as they exhibit entirely different error sources. They include pathintegral Quantum Monte Carlo (PIQMC) calculations, 3,4 Density Matrix Quamtum Monte Carlo (DMQMC) calculations, 44 The Kohn−Luttinger conundrum is also closely related to this work. It states that the infinite-size zero-temperature limit of finite-temperature many-body perturbation theory not necessarily agrees with the infinite-size limit of zero-temperature many-body perturbation theory. In the common approach where the zero-momentum part of the electrostatic interaction is disregarded, certain terms called anomalous diagrams affect both, the chemical potential and the grand potential. It has been shown that their contributions cancel in the zero-temperature limit of the free energy under isotropic conditions 7 and later more generally. 53,54 Hirata and co-workers developed a perturbation theory that simultaneously produces order-byorder corrections to the grand potential and the chemical potential, while fixing the expected number of electrons. 55,56 Employing the number-conserving perturbation theory and treating possible degeneracies of the SCF reference in the lowtemperature limit, 57 they overcome the conundrum. 58 With the method of this work the situation is different and will be discussed at the end of subsection 2.4. Further discussions can be found in refs 18−20 and 59.

METHODS
Let us now develop the regularization approach for the prototypical warm-dense system: the warm uniform electron gas (UEG). The UEG is a model of a metal, where the positive ions of the lattice are replaced by a static homogeneous positive background charge. It has a vanishing band gap in the infinitesize limit and thus qualifies for a warm-dense system at all nonzero temperatures. All properties of the warm UEG depend only on the thermodynamic state, specified by its density and temperature. The density is usually given in terms of the Wigner−Seitz radius r s in atomic units, such that the volume of a sphere with radius r s corresponds to the average volume per electron. It is also convenient to specify the temperature in terms of the dimensionless ratio θ = k B T/ε F , where k B T is the average thermal energy and ε F = k F 2 /2 is the Fermi energy of a free spinunpolarized (paramagnetic), infinite electron gas at the corresponding density and at zero temperature with k F 3 = 9π/ 4r s 3 . This defines a natural temperature scale where different densities can be compared to each other more directly.
The UEG is modeled by a finite cubic box of length under periodic boundary conditions having the volume = = r 4 /3 3 s

3
. It contains a homogeneous positive charge density with a total charge of elementary charges, which is considered fixed. In the grand-canonical ensemble the number of electrons in the system is not fixed but rather fluctuates around its expectation value which depends on the chemical potential μ. Later, μ will be chosen such that the expected number of electrons N N equals, or is close to, the number of positive charges . To treat the diverging electrostatic interaction, the Spencer−Alavi or the Yukawa regularization of the electrostatic interaction is used. The Spencer−Alavi interaction is given by the usual Coulomb interaction 1/r 12 for the distance between two electronic coordinates < r 12 that lie within a sphere of radius . The interaction is zero otherwise. The truncation radius = r s 1/3 is coupled to the system size and the simultaneous limit 3 is studied. The kernel of this interaction within a sum over momenta is = V q q q ( ) 4 (1 cos )/ 2 . For finite the kernel is also finite at q = 0 and evaluates to 2 / 2 . With this choice the interaction "sees" on average electrons. In the limit it reduces to the usual electrostatic interaction 1/r 12 for all distances. 9 The Yukawa regularization is also studied for comparison. It is given by exp(−αr 12 )/r 12 in real space and by = + V q q ( ) 4 / ( ) 2 2 in momentum space in a sum over states. The regularization parameter α is chosen such that the weighting function exp(−αr 12 ) integrates to , which means that each electron "sees" electrons. With this choice, the regularization parameter is coupled to the system size satisfying = r 6/ 3 s 3 and the limit 3 is studied. For numerical computations the Yukawa regularization is an inefficient choice since asymptotic behavior usually only sets in at impracticably large system sizes. While it is not useful for retrieving absolute free energies, the difference between free energies reaches asymptotic behavior for the system sizes studied in this work. Figure 1 illustrates the regularized electrostatic interaction energy between unit charges at a distance r 12 for the employed regularization schemes for different system sizes.
All states are expanded in antisymmetrized products of oneelectron wave functions that are eigenfunctions of the singleelectron kinetic operator −∇ 2 /2 under periodic boundary conditions. The normalized eigenfunctions are the plane waves commensurate with the box length It consists of three terms: the kinetic term T, the electron− background interaction V ext and the electron−electron interaction V , respectively. Exact diagonalization of the Hamiltonian is infeasible except for very limited system sizes. This work shall employ the approximation approach of computational materials science, where one first performs a self-consistent field calculation, followed by an approximation of the correlation based on the SCF result. In accordance with the common workflow of random phase approximation (RPA) calculations for low-band gap systems, the SCF only employs the Hartree approximation rather than Hartree−Fock and exchange is considered at first order non-self-consistently. 60 The finite temperature correlation contributions are estimated by a linearized form of the direct-ring coupled-cluster doubles approximation.

Self-Consistent Field in the Hartree Approximation.
In the self-consistent field approach the two-body operator in the electron−electron interaction V is partially contracted to a one-body interaction. 37 In the Hartree approximation only the direct contraction is considered and the resulting one-body operator is given by where A 0 denotes the one-body thermal equilibrium expectation value of the operator A, defined by with the (non-normalized) one-body density matrix Figure 1. Illustration of the Spencer−Alavi (red) and the Yukawa (orange) regularization of the electrostatic interaction for different system sizes . The respective regularization parameter or α is coupled to the system size such that the scope of the interaction agrees with the system volume .
All terms in eq 3 are diagonal in the chosen basis so we can immediately write the equations for the eigenvalue of each state i = (k i , σ i ) Introducing the notation ε i = k i 2 /2 + Δε, we have to find a shift of eigenenergies Δε, uniform for all states, satisfying the nonlinear equation for the given thermodynamic state point ( , , ). So far, the number of positive charges is an independent parameter. The quantity ΔN 0 denotes the net-negative charge of the system in the noninteracting approximation. Note that it may differ from zero for charge neutral systems as the fully interacting = N N may differ from its noninteracting approximation = N N 0 0 . Having solved eq 7 for Δε we can evaluate the noninteracting grand potential with the ef fective chemical potential η ≔ μ − Δε. If we want to compare energies per electron for different system sizes we also need to account for the background−background interaction energy, which is independent of the electronic degrees of freedom. Furthermore, pairwise interactions in Ω 0 are doublecounted in the SCF. Accounting for both contributions yields the mean-field grand potential in the Hartree approximation Note that the expected number of electrons in the Hartree approximation N H = −∂ μ Ω H equals the expected number of electrons of the noninteracting system N 0 = −∂ μ Ω 0 , which is given by the sum ∑ i n i of noninteracting occupancies The Hartree grand potential Ω H (μ) is a function of the chemical potential μ. Of particular interest is the chemical potential μ H for which the expected number of electrons in the Hartree approximation N H matches the number of positive charges . This chemical potential is chosen for determining the Hartree free energy from a Legendre transformation For this particular chemical potential, the solution of the Hartree equation Δε = 0 follows trivially from eq 7 and the effective chemical potential η H = μ H − Δε equals the Hartree chemical potential.
A Rough Estimate of the Hartree Self-Consistent Field. Before turning to the other contributions to the grand potential, it is interesting to estimate the Hartree solution Δε for chemical potentials μ close to the chemical potential μ H , which satisfies = N H . To this end, we expand the expected number of electrons N H = −∂ μ Ω 0 as a function of the effective chemical potential η = μ − Δε at the Hartree effective chemical potential between the expected number of electrons in the SCF calculation and the number of positive charges to first order in (η − μ H ) by We further assume that in the warm-dense-matter regime the term ∑ i n i (1 − n i ) scales linearly with system size. Together with eq 7 Δε = V(0)ΔN 0 we are now in the position to approximately solve eq 11 for ΔN 0 : Inserting ΔN 0 into eq 7 and expanding in powers of for large finally yields for the Spencer−Alavi truncated Coulomb kernel. Interestingly, the ( ) 1/3 terms in the above equations are exactly linear in (μ − μ H ), although ΔN 0 may have terms depending on higher orders of (μ − μ H ). This stems from the relation N 1/3 0 , which is exactly linear. In case of the Yukawa regularization, V(0) is 4 / 2 which gives a prefactor of 6 2/3 r s /3 instead of 2r s /3 in both of the above equations.
This is an essential result. For a finite deviation of the chemical potential μ from the noninteracting chemical potential μ H = μ 0 , the expected number of excess electrons N 0 of the entire system diverges as ( ) 1/3 . It grows proportionally to the range of the regularized Coulomb interaction. However, the expected number of excess electrons per positive charge always converges to one as ( ) 2/3 and the electrostatic energy of the system per electron Δε asymptotically reflects the change of the chemical potential. When changing the chemical potential, the relaxed self-consistent field eigenenergies asymptotically change in the exact same way for large enough system sizes. The numerical results in subsection 3.1 indicate that this behavior already sets in for relatively low system sizes. Since only the difference of the chemical potential μ and the eigenenergies ε i = k i 2 /2 + Δε enter in subsequent calculations, size-intensive observables will be asymptotically independent of the choice of μ in the thermodynamic limit. At warm-dense conditions, any finite choice is acceptable, including μ = 0, which is the classical definition of the chemical potential of electrons in a grounded conductor that can supply or absorb any number of electrons.
Decoupling the Infinite-Range Limit of the Coulomb Interaction. Instead of the simultaneous limit of 3 , studied in this work, one could also perform the infinite-range limit before performing the infinite- respectively. With increasing range , deviations become increasingly expensive and one arrives asymptotically at the canonical ensemble for exactly = N 0 electrons, where the chemical potential μ plays no role anymore. However, classical systems of charged particles exhibit larger finite-size errors when simulating in the canonical ensemble compared to simulating in the grand-canonical ensemble with the range of the Coulomb interactions coupled to the size of the simulation cell. 8 Thus, this work pursues the latter approach.
Using Fixed Instead of Relaxed Orbitals. If one chooses to work with a fixed set of orbitals and eigenenergies for various values of the chemical potential μ, the electron−background interaction and the electron−electron density interaction do not cancel for μ deviating from μ H . There are terms of order n in the perturbation series whose contributions scale as ( ( ) ) n 2/3 . They are not extensive for each n and alternating in sign. This can be cured by summing the interactions to infinite order, which is equivalent to performing an SCF calculation at the modified chemical potential. 61 Fluctuations This agrees qualitatively with classical charge fluctuations δN C 2 on the surface of a grounded conducting sphere of radius = r s 1/3 , found from the equipartition theorem Note that we need to consider the response of the one-body energies to changes of the chemical potential for computing derivatives of the grand potential beyond first order, such as ∂N H /∂μ = −∂ 2 Ω 0 /∂μ 2 in eq 15. For comparison, a fixed density matrix 0 that separates into a product of one-body density matrices is only capable of describing electron-number fluctuations of the form , which is proportional to rather than to 1/3 at warm-dense conditions.

First-Order Exchange.
The self-consistent field approximation is crude but computationally efficient. To improve on the approximation, finite temperature many-body perturbation theory (FT-MBPT) offers an expansion of the grand potential in powers of the difference between the true Hamiltonian Ĥand the self-consistent field Hamiltonian H 0 . Having employed the Hartree approximation for the SCF, the leading order term is the first-order exchange term where we again use compound indices i = (k i , σ i ) and j = (k j , σ j ) to denote the spatial and spin components of the respective spinorbitals. We also employ the shorthand notation n ij... = n i n j ... for products of one-body occupancies. V sr pq denotes the components of the electron−electron interaction operator in the basis of the plane-wave spin-orbitals such that For a translationally invariant, isotropic interaction the components read The first-order exchange term with non-Hartree−Fock orbitals is often referred to as exact exchange (EE). It is given by (19) Adding the first-order exchange contribution Ω x to Ω H yields the improved Hartree-exchange approximation Ω Hx .

Linearized Direct-Ring Coupled Cluster.
Let us now turn to correlation and exchange effects beyond first order. Here, it is treated at the level of linearized direct-ring coupled-cluster doubles (ldrCCD) theory. 30 A truncation of the perturbation expansion at any finite order n diverges as T ( ) n ( 1) for the uniform electron gas for low temperatures T in the infinite-size limit. On the other hand, summing over the so-called ring terms up to infinite order n yields convergent results at zero temperature in the limit after n → ∞. 62,63 The analogous ring-term summation of finite-temperature MBPT also yields convergent results for all temperatures although its limit T → 0 might differ from the result obtained from a zerotemperature theory. We desire a theory with such a convergence behavior for T → 0, at least in principle. ldrCCD is one of the simplest theories providing this resummation of the ring terms. It contains all ring terms that can be formed with exactly two particle/hole pairs and additionally contains their corresponding screened-exchange terms. It is determined by the finitetemperature linearized direct-ring coupled-cluster amplitude integral equations  (20) with Δ ij ab = ε a − ε i + ε b − ε j and where we now also need products of vacancy and occupancy probabilities, denoted by = n n n (1 ) with = n n n n n (1 ) (1 ) ij ab a i b j . All indices iterate in principle over the infinite number of plane wave states. Practical truncation schemes are discussed in section 3.3. The linear system of coupled integral equations in eq 20 can be solved by diagonalizing an effective particle/hole interaction H̃, analogous to the Tamm−Dancoff approximation of the Casida equations at zero temperature. The effective particle/hole interaction reads (22) which, interpreting the indices (b, j) as a compound row index and the indices (a, i) as a compound column index, is a hermitian matrix and thus permits a real-valued eigendecomposition. We can then transform the electron repulsion integrals with and without exchange into the space of eigenmodes  (24) and finally retrieve the ldrCCD approximation of the correlation grand potential from i k j j j j j y with Λ FG = Λ F F + Λ G G and * = V V FG FG denoting the conjugate transpose. 30 2.4. Free Energies. So far, we have discussed all considered contributions to the grand potential as a function of the thermodynamic state point in the grandcanonical ensemble, in particular of the chemical potential μ. The number of positive charges is merely a system parameter. We are, however, interested in the free energy F ( ) Hxc of the charge-neutral system where the expected number of electrons N Hxc ≔ −∂ μ Ω Hxc equals the fixed number of positive charges. It is found from the Legendre transformation Hxc Hxc Hxc (27) where μ Hxc satisfies the charge-neutrality condition for the H a r t r e e -e x c h a n g e -c o r r e l a t i o n g r a n d p o t e n t i a l = ( ) to eq 14. From eqs 19 and 25 it follows that the only terms that depend on μ in the remaining contribution Ω xc are the o c c u p a n c y a n d v a c a n c y e x p e c t a t i o n v a l u e s and n a = 1 − n a , respectively. The expectation values n i depend only on the difference ε i − μ between the eigenenergies ε i = k i 2 /2 + Δε and the chemical potential μ, where Δε is the shift of eigenenergies, uniform for all states i, found from the self-consistent field solution for the interacting chemical potential μ Hxc . Using the notion of the effective chemical potential η = μ − Δε introduced in subsection 2.1, we can write the derivative with respect to the chemical potential in terms of a derivative with respect to the effective chemical potential from the chain rule xc xc (28) We have already estimated the asymptotic behavior of Δε in eq 13 from which we can find the behavior of ∂ μ η(μ) for large : Note that −∂ ημ 2 Ω 0 = β∑ i n i i scales linearly with under warmdense conditions. Similarly, since ∂ η n a = −βn a a , we can also assume that −∂ η Ω xc scales at most linearly with the system size under these conditions. Collecting all contributions to the expected number of electrons gives Although the expected number of electrons per positive charge converges to one for any chemical potential in the thermodynamic limit, there is a nonvanishing deviation from the noninteracting chemical potential μ H required if also the absolute expected number of electrons N Hxc should match the number of positive charges for large . Knowing the asymptotic behavior of the interacting chemical potential μ Hxc we can now estimate the free energy for large system sizes. For that purpose, we expand the Hartree-exchange-correlation free energy at the noninteracting chemical potential in eq 27 in powers of the difference (μ Hxc − μ H ), which we have found to be finite but approximately independent of : Our estimate of ∂ μ η in eq 29 is approximately independent of μ. Therefore, the higher derivatives of the exchange-correlation grand potential occurring in the above expansion are estimated to be of the form ( )( ) In the thermodynamic limit the exchange-correlation grand potential per electron, evaluated at the noninteracting chemical potential, is estimated to agree with the exchange-correlation free energy per electron = f F lim ( )/ xc xc , found at the interacting chemical potential. As for the SCF estimates, using the Yukawa regularization yields a prefactor of 6 2/3 r s /3 instead of 2r s /3 in the above equation.
This is the main result of this work, and the numerical studies in the following section show that this asymptotic estimate applies already at relatively small system sizes in the uniform electron gas for the densities and temperatures considered. For finite system sizes , eq 34 relates the difference between Ω xc (μ H ) and F ( ) xc to the difference between the interacting and the noninteracting chemical potential. The latter converges faster with system size and this relation permits an estimate of the remaining finite-size error in Ω c for the thermodynamic limit extrapolation. One can also employ eq 31 to estimate the interacting chemical potential μ Hxc from a correlation calculation at a noninteracting chemical potential if the derivative with respect to η can be found efficiently.
Comparison with the Madelung Technique and the Kohn−Luttinger Conundrum. This work examines the simultaneous limit of the interaction range and the system size 3 . The scope of the electrostatic interaction is limited to the size of the system . This is the key ingredient for the free energy per electron to be asymptotically independent of the chemical potential. A spherically truncated interaction and the Yukawa regularized interaction are studied in detail in this work. Other truncation schemes, such as the minimal image convention or the Wigner−Seitz truncation, are also expected to yield this behavior, albeit with different prefactors. If the infiniterange limit is taken before the thermodynamic limit , the ensemble reduces to the canonical ensemble, as discussed at the end of subsection 2.1, which is expected to exhibit larger finite-size errors already in classical systems. 8 A common alternative to limiting the scope of the electrostatic interaction is the Madelung technique. It considers the interaction energy ξ < 0 of an elementary point charge with its own periodic images and a neutralizing background charge for the zero-momentum part of the Coulomb interaction V(q = 0) = −ξ. For other momenta q with q = |q| > 0, the standard Coulomb interaction = V q q ( ) 4 / 2 is employed. 64 This technique was developed for Monte Carlo simulations in the canonical ensemble and assumes a charge-neutral simulation cell under periodic boundary conditions. However, it has also been applied to finite temperature MBPT calculations in the grand-canonical ensemble, 28,49 where charge-neutrality occurs only on average unless the fluctuations of electrons and positive charges are fully correlated�an unlikely scenario considering their different masses. Otherwise, this choice for the Hamiltonian is an approximation. In the limit of infinitely large simulation cells, the Madelung term vanishes ξ → 0 and one arrives at the assumption V(q = 0) = 0 that Kohn and Luttinger make to arrive at eq 20 in ref 7. In this work they find that the infinite-size and zerotemperature limit of the grand potential in metals depends on the order of the limits and reduces to the value produced by zero-temperature MBPT only for after T → 0 and if the zero-temperature limit of the SCF reference is nondegenerate. In the paramagnetic case, all orbitals of a common eigenenergy must be asymptotically doubly occupied or not at all. For the simple-cubic UEG this only applies to system sizes of = 2, 14, 38, 54, ... This discontinuity is called Kohn− Luttinger conundrum, although the authors never used this term. In the reverse order, the result differs by a non-zero value, stemming from so-called anomalous terms. 7 Such terms contain orbitals that are simultaneously particles and holes and these terms are nonvanishing if the density of states at the Fermi surface is nonvanishing. For an electron−electron interaction satisfying the assumption V(q = 0) = 0 they also show that the correlation contribution of the anomalous terms to the chemical potential exactly counteracts the anomalous terms in the grand potential when regarding the Legendre transformed free energy for isotropic systems.
With the method proposed in this work, the situation is different. For a fixed finite T > 0, the free energy per electron is asymptotically independent of the chemical potential for according to eq 34. Correlation corrections to the chemical potential play no role for the free energy per electron and contributions from anomalous terms to the exchangecorrelation grand potential are not canceled. Taking the limit T → 0 after may yield a different result than after T → 0 in metallic systems.
In this work, only finite systems have been considered that have no partially occupied orbitals in the limit T → 0. For each system size, the conventional closed-shell zero-temperature ldrCCD calculation converges and agrees with the T → 0 limit of the finite-temperature ldrCCD calculation at the respective system size. 30 If the finite-temperature SCF reference becomes degenerate for T → 0, as in system sizes deviating from the simple-cubic closed-shell numbers = 2, 14, 38, 54, ..., finiteorder thermal MBPT calculations diverge, as demonstrated by Hirata. 57 ldrCCD calculations have been shown to converge to a finite value for such systems owing to the infinite-order resummation. 30 Still, it remains to be studied whether an infinite-size extrapolation using degenerate system sizes, which deviate from the above closed-shell numbers, yields the same Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article result as the thermodynamic limit extrapolations from closedshell system sizes presented in this work. At zero temperature, it is a common practice only to choose closed-shell sizes.
Extrapolations of finite-system calculations to the infinite-size limit agree well with continuous momentum integration-based methods that are in the thermodynamic limit by construction. 65,66 At finite-temperature, only few resummed theories, such as the random phase approximation (RPA), 24 can be integrated in momentum-space in practice, allowing the study of the T → 0 limit after the limit in metallic systems. However, the only anomalous terms contained in direct-ring theories, such as RPA and ldrCCD, are particle/hole rings with zero momentum transfer q = 0. Including or excluding this point has no effect on the integral. Yet, there are more anomalous terms and whether their sum in other theories is zero, finite, or even diverging in the grand potential is an open question and topic of ongoing research. An alternative question is, how the number-conserving perturbation theory of Hirata and co-workers 55,56,67 can be integrated in an infinite-order theory, such as RPA or coupledcluster.

NUMERICAL RESULTS
To assess the large system-size estimates in the previous section, numerical calculations of the paramagnetic uniform electron gas have been conducted for simple-cubic systems containing 38,54,66,114,162,246,294,342,358, and 406 electrons on average. The system sizes have been chosen such that degenerate spatial orbitals can be fully occupied at zerotemperature in a closed-shell self-consistent field calculation.
3.1. Hartree Self-Consistent Field. Unlike in Hartree− Fock, where exchange is included in the SCF calculation, the metallic spectrum of the free uniform electron gas remains in the Hartree SCF. No band gap forms since the eigenenergies of mostly occupied and mostly unoccupied orbitals are affected in the same way. At zero temperature this has been found to be relevant for RPA-like calculations. 60 Consequently, each eigenvalue in eq 6 only depends on its kinetic energy and the sum of all occupancies. A uniform shift of the eigenenergies Δε is the only number that needs to be found, although in a nonlinear equation. At finite temperature, all orbitals contribute in principle. In this work the number of spatial orbitals for the closed-shell SCF calculation has been truncated at roughly 800 times the number of orbitals occupied at zero temperature. Sums over the orbitals beyond this number occurring in Ω 0 and N H have been approximated by integrals. With this treatment, all SCF quantities are well converged and the computation time for the SCF calculation is still negligible compared to the correlation calculations. Note that the SCF calculations have been repeated to yield relaxed eigenenergies for each value of the chemical potential in search for the chemical potential μ Hxc where the expected number of electrons matches the number of positive charges .
Only the difference between the eigenenergies ε i = k i 2 /2 + Δε and the chemical potential μ occurs in the expressions of manybody perturbation theory where Δε depends on μ. Thus, they can be viewed rather as functions of the effective chemical potential η = μ − Δε. Figure 2 shows how the effective potential changes when the chemical potential is changed. It plots the derivative ∂η/∂μ against 2/3 where is the system size. The derivative has been evaluated at the fully interacting chemical potential μ Hxc , except for the largest system size = 23674, where no correlation calculation has been conducted and μ H has been used instead. Already for moderate system sizes, a change of the chemical potential has about 2 orders of magnitude less an effect on η and in consequence on the expressions of FT-MBPT. For large the effect on η decreases, scaling as ( ) 2/3 , as estimated in eq 29, and vanishes in the thermodynamic limit.

First-Order Exchange.
The exchange contributions to the grand potential Ω x have been evaluated according to eq 19 using all orbitals that have been considered in the SCF calculation. The convergence with the number of closed-shell orbitals is faster than for the SCF quantities and no analytic treatment of the orbitals beyond 800 times the zero-temperature orbitals is necessary. Although considerably more demanding computationally than the SCF calculation, its evaluation is still negligible compared to the correlation calculation. The derivative of the exchange contribution with respect to η for the expected number of electrons has been evaluated analytically.
3.3. Linearized Direct-Ring Coupled Cluster. The correlation and exchange effects beyond first order have been approximated on the level of closed-shell linearized direct-ring coupled cluster doubles (ldrCCD) theory. It is one of the simplest theories whose zero-temperature and infinite-size limit exists. Still, it is expected to capture the dominant part of the long-range correlation. The advantage of ldrCCD is that it can be evaluated from a diagonalization of an effective particle/hole Figure 2. The effective potential η = μ − Δε is a measure for how quantities in the perturbation expansion depend on μ. The figure shows that the change of the effective potential η with respect to the chemical potential μ decreases with increasing system size. An extrapolation of the largest calculations with 2/3 indicates that the effective potential η becomes independent of μ in the thermodynamic limit.
Hamiltonian and consequently permits an analytic imaginary time integration. Apart from numerical considerations, the temperature can be arbitrarily low. Let N p denote the number of spatial orbitals considered for the ldrCCD calculation. The direct-ring structure of the effective Hamiltonian in eq 22 is momentum conserving. In a uniform system and due to point-group symmetry it is therefore sufficient to consider independent N p × N p matrices for each momentum difference vector q = k b − k j in the wedge 0 ≤ q x ≤ q y ≤ q z instead of one N p 2 × N p 2 matrix. For the largest system size 114 independent 4385 × 4385 matrices have been diagonalized. The matrices are real valued and symmetric and can be diagonalized efficiently with standard linear algebra packages.
Unlike at zero-temperature, the spectrum of each matrix is not necessarily positive-definite. Negative eigenvalues Λ F F can occur when eigenenergies of contributing hole-orbitals are above the energies of contributing particle-orbitals, which is possible at finite temperature. Negative eigenvalues pose numerical difficulties occurring in the exponent of eq 25. However, the final product with the square roots of the occupancy products n ij ab in eqs 23 and 24 leads to a finite contribution. In practice, the term δ a b δ j i Δ j b has been truncated to zero if the occupancy and vacancy product n j b was below 10 −12 . For the ldrCCD calculation N p has been chosen about 20 times the number of zero-temperature occupied spatial orbitals. The contribution from the orbitals beyond has been extrapolated from the asymptotic behavior of RPA-like correlation energies. The finite-basis-set error scales as q ( ) max 3 , where q max is the magnitude of the largest considered plane wave momentum difference. 68 A Hann window has been used to obtain a soft cutoff for four different values of q max to smooth the samples for the q max −3 extrapolation to the complete basis set (CBS) limit. 60 The correlation coefficients of the regression curves range between 0.97 and practically 1. The 67% confidence intervals of the CBS limits are given in the ±CBS column in Table 1.
For each system size, multiple calculations of Ω Hxc have been conducted in search for the chemical potential μ Hxc where N Hxc = −∂ μ Ω Hxc agrees with the number of positive charges . The derivative of Ω c has been evaluated numerically from a polynomial fit. The next estimate Hxc at the current chemical potential μ has been found from the difference of N Hxc assuming that the dominant change in N Hxc stems from the change in N H = −∂ μ Ω 0 . This gives an equation for the dominant change in the chemical potential where all involved quantities can be readily evaluated at the current chemical potential μ. This procedure has required about 8 iterations until convergence for each considered system size . All energies are given in Hartree. f xc is retrieved from the thermodynamic limit and complete-basis-set limit of Ω xc . Exchange and correlation contributions have been extrapolated separately. The expected statistical errors from the infinite-size and infinite-basis-set extrapolations of the correlation contributions are given in the ±TDL and ±CBS column, respectively. The TDL and CBS errors of the exchange contributions are negligible. For comparison, the free energies obtained from the improved Singwi−Tosi−Land−Sjolander 70 parametrization of Quantum Monte Carlo results and the difference to ldrCCD are given in the columns f xc iSTLS and Δf xc , respectively. The left panel of Figure 3 depicts how the various contributions to the exchange-correlation grand potential Ω xc depend on the reduced temperature for a fixed density of r s = 2 and two system sizes = 162 and = 246. The right panel shows the temperature dependence of the chemical potentials satisfying the charge-neutrality condition at increasingly accurate levels of theory. The = 162 and = 246 results are denoted by pluses (+) and crosses (×), respectively. Deviations between these system sizes at low temperatures originate from a partition of the orbitals into (almost) purely occupied or purely unoccupied orbitals, occurring at lower temperatures for larger systems. A T → 0 extrapolation has not been conducted.
3.4. Free Energies in the Thermodynamic Limit. Finding the thermodynamic limit poses a difficult task in the calculation of extended systems. First, we assess whether the asymptotic behavior estimated by eq 34 applies in the UEG as a prototypical warm-dense system. Figure 4 plots the difference between the exchange-correlation free energy per electron and the exchange-correlation grand potential per electron, evaluated at the Hartree chemical potential μ H , against the system size 2/3 . In this graph, an asymptotic behavior as estimated from eq 34 is expected to appear as a line through the origin. As a guide to the eye, the results are connected with dashed red lines. The linear extrapolations from the largest system sizes are shown as solid red lines. The 67%-confidence intervals of the thermodynamic limits are indicated by the error bars on the vertical axis. They confirm numerically that the two exchangecorrelation free energies agree in the thermodynamic limit of the warm UEG for all densities and temperature considered.
The terms in Ω xc converge with different rates to the thermodynamic limit. At the largest considered system sizes the exchange contributions are almost converged. The remaining correlation terms converge as ( ) 2/3 in the low temperature . Finite-size dependence of the difference between the exchange-correlation free energy per electron and the exchange-correlation grand potential per electron for various densities and temperatures. The correlation contributions are approximated by the linearized direct-ring coupled cluster doubles (ldrCCD) theory. The grand potential has been evaluated at the noninteracting chemical potential μ H while the free energy requires the correlated chemical potential μ Hxc . Extrapolations of the largest system sizes show that the two free energies coincide in the infinite-size limit. regime and as ( ) 1 otherwise. 51 Also, the effective chemical potential η, which Ω c depends on, converges as ( ) 2/3 . Thus, F xc and Ω xc are also individually expected to converge to the thermodynamic as ( ) 2/3 . Figure 5 plots F xc and Ω xc individually against the system size 2/3 , while Figure 4 showed only their difference. Unlike their difference, the two energies suffer from shell effects. They could be alleviated by twist averaging 51,69 but this has not been done in this work. The solid lines show the 2/3 fit for the largest system sizes of the respective sets and the statistical error of the infinite-size extrapolation for both energies is indicated by the error bars on the vertical axes. Interestingly, in most cases the slope of the grand potential extrapolation is flatter than that of the free energy extrapolation, making the extrapolation of the grand potential less dependent on the functional form of the asymptotic behavior. Table 1 summarizes the exchange-correlation free energies f xc found in the thermodynamic limit and gives the 67% confidence interval of the infinite-size extrapolation in the ± TDL column. At low temperatures, the results compare well to the improved Singwi−Tosi−Land−Sjolander parametrization 70 of previous Quantum Monte Carlo calculations despite the simple ldrCCD theory employed. For zero-temperature this fortuitous agreement of ldrCCD and also of drCCD has already been shown by Freeman. 65 Reference 51 gives an overview over other parametrizations. At higher temperatures, however, the difference is sizable and ldrCCD underestimates the magnitude of the negative exchange-correlation free energy considerably. Figure 5. Finite-size dependence of the exchange-correlation grand potential Ω xc per electron in red and the exchange-correlation free energy F xc per electron in green for various densities and temperatures of the nonspin-polarized (paramagnetic) warm dense electron gas. The difference of the two graphs does not exhibit shell effects and is shown in Figure 4. The finite-size error of the exchange contribution is negligible for the largest system sizes and the remaining terms are expected to converge as ( ) 2/3 to the thermodynamic limit. The solid red and green lines show the extrapolations to the thermodynamic limit of the grand potential and the free energy, respectively. The error bars on the vertical axes indicate the statistical error of the extrapolations.
All of the previous numerical results have been produced employing the Spencer−Alavi spherical truncation scheme for all electrostatic interactions: the background−background interaction, the electron−background and the electron− electron interaction in the Hartree SCF, as well as the electron−electron interaction in the perturbation. To assess to what degree the estimates of the asymptotic behavior depend on the employed regularization scheme, the Yukawa regularization has been implemented as well and calculations have been conducted for various densities and temperatures. It can bee seen in the lower two panels of Figure 6 that the exchangecorrelation free energy and the exchange-correlation grand potential individually hardly converge to the asymptotic domain using the Yukawa regularization. Note the different energy scales. Although only small shell effects are present, as opposed to the spherical truncation scheme, the Yukawa regularization is not useful in practice for retrieving converged results in the thermodynamic limit. In contrast, the difference of the two energies already exhibits an asymptotic behavior that agrees with the estimated behavior.

SUMMARY
This work shows that the infinite-size limit of finite temperature many-body perturbation theory can be found efficiently with a truncated Coulomb interaction. The truncation radius is chosen such that the volume of the interaction agrees with the volume of the simulated cell. Such schemes have previously been studied in classical systems of electrostatically interacting particles, as well as for the Fock-exchange contribution in zero-temperature MBPT. Here, the truncation scheme is employed for all electrostatic interactions in the uniform electron gas: electron−electron, electron−background, and background− background.
It is found that due to the long-ranged nature of the electrostatic interaction the difference of the average number of mobile electrons and fixed positive charges scales asymptotically as ( ) ( ) Hxc 1/3 for large system sizes where is the number of positive charges and μ Hxc is the chemical potential where the number of electrons equals the number of positive charges, including exchange and correlation effects. Thus, the ratio of the number of electrons and positive charges tends to one for any finite choice of the chemical potential μ.
An important consequence is that also the exchangecorrelation grand potential per electron, evaluated at the noninteracting Hartree self-consistent field chemical potential μ H , asymptotically agrees with the free energy per electron, found from a Legendre transformation at the interacting chemical potential μ Hxc : The latter requires multiple iterations of the expensive correlation calculations during the nonlinear search for the interacting chemical potential for each system size considered for the thermodynamic limit extrapolation. The above asymptotic behavior has been estimated in general for matter under warm-dense conditions and it has been shown explicitly for the warm uniform electron gas for various densities Figure 6. Top panels show the finite-size dependence of the difference between the exchange-correlation free energy and the exchange-correlation grand potential for θ = 0.5 and two densities employing the Yukawa regularization, rather than the Spencer−Alavi truncation. The bottom panels show these free energies individually. Although neither the exchange-correlation free energy nor the exchange-correlation grand potential have reached the asymptotic domain individually in case of the Yukawa regularization, their difference has�just barely: Only the largest 3 system sizes have been used for the extrapolation. and temperatures employing the linearized direct-ring coupled cluster doubles theory for approximating exchange-correlation effects. The author declares no competing financial interest. The program and the settings used to produce the data in this work are publicly available at https://gitlab.cc4s.org/cqt/wegldrccd.git

■ ACKNOWLEDGMENTS
The author wishes to thank Isabella Floss, Andreas Irmler, Evgeny Moermann, Nikolaos Masios, Andreas Savin, Sam Trickey, and Corbinian Wellenhofer for constructive discussions and remarks on the manuscript. The author gratefully acknowledges TU Wien Bibliothek for financial support through its Open Access Funding Programme and the group of Andreas Gruneis for providing cpu time on its computational resources.