How the Exchange Energy Can Affect the Power Laws Used to Extrapolate the Coupled Cluster Correlation Energy to the Thermodynamic Limit

Finite size error is commonly removed from coupled cluster theory calculations by N–1 extrapolations over correlation energy calculations of different system sizes (N), where the N–1 scaling comes from the total energy rather than the correlation energy. However, previous studies in the quantum Monte Carlo community suggest an exchange-energy-like power law of N–2/3 should also be present in the correlation energy when using the conventional Coulomb interaction. The rationale for this is that the total energy goes as N–1 and the exchange energy goes as N–2/3; thus, the correlation energy should be a combination of these two power laws. Further, in coupled cluster theory, these power laws are related to the low G scaling of the transition structure factor, S(G), which is a property of the coupled cluster wave function calculated from the amplitudes. We show here that data from coupled cluster doubles calculations on the uniform electron gas fit a function with a low G behavior of S(G) ∼ G. The prefactor for this linear term is derived from the exchange energy to be consistent with an N–2/3 power law at large N. Incorporating the exchange structure factor into the transition structure factor results in a combined structure factor of S(G) ∼ G2, consistent with an N–1 scaling of the exchange-correlation energy. We then look for the presence of an N–2/3 power law in the energy. To do so, we first develop a plane-wave cutoff scheme with less noise than the traditional basis set used for the uniform electron gas. Then, we collect data from a wide range of electron numbers and densities to systematically test five methods using N–1 scaling, N–2/3 scaling, or combinations of both scaling behaviors. We find that power laws that incorporate both N–1 and N–2/3 scaling perform better than either alone, especially when the prefactor for N–2/3 scaling can be found from exchange energy calculations.


INTRODUCTION
There has been a recent push toward developing wave functionbased methods such as coupled cluster theory for solids. 1−21 A long-term goal of this work is to provide highly accurate energy calculations for materials design. Coupled cluster has been growing in popularity for solid state calculations due to its ability to obtain the correlation energy (i.e., E total − E HF , where E HF is the Hartree−Fock energy) in a versatile and systemically improvable way. However, one of the main issues facing coupled cluster is that the energies show slow, polynomial scaling as a function of both system size, N, and k-points when converging to the thermodynamic limit (TDL)�the limit of an infinite atom or particle number. As most energy calculations gain meaningful insight about the system at the thermodynamic limit, it is imperative that we know the exact rate at which the coupled cluster correlation energies approach this limit.
Recent advances in coupled cluster theory have made coupled cluster single and doubles (CCSD) calculations for solids seem increasingly routine, 17−21 overcoming numerical convergence issues with small denominators, the divergences of perturbative methods, and technological barriers. In our previous studies, we have found that CCSD is a reliable way to study finite size effects (FSE) for the coupled cluster hierarchy of methods, especially when basis set errors can be effectively controlled. 17,22−26 In turn, the study of finite size effects is important in ensuring that coupled cluster theory is generally useful for energy calculations of solids.
A popular way to address the cost scaling issue and obtain TDL energy estimates from coupled cluster calculations for smaller system sizes is to perform an extrapolation to the TDL. The TDL-extrapolated energy is typically calculated by running increasingly large system sizes and then fitting the energies at the larger N to the following function: E N = lim N→∞ (E TDL + mN −γ ).
Here, N is the system size and refers to a number of electrons. The number of k-points, N k , can also be used. The variable γ defines the convergence rate. If γ is known exactly for large N, the TDL energy can be estimated more accurately. For other energies, such as the correlation energy, there is an ongoing discussion in the literature as to the exact value of both γ and the form of the extrapolation equation itself.
A commonly assumed convergence rate for the correlation energy is N −1 , 7,19,27 the same as the total energy relationship.
The N −1 convergence of the correlation energy, which is physically attributed to long-range van der Waals forces 7,28,29 and can be derived in the UEG, 30 has been related to the low momentum limit (G → 0) of the transition structure factor S(G). The transition structure factor is calculated from the amplitudes of the CCSD wave function, and the sum over its pointwise product with the Coulomb operator in k-space yields the correlation energy. As such, its scaling at low momenta relates to the power law of the TDL extrapolation: a convergence of S(G) ∼ G 2 predicts a power law of N −1 . 23,27,31−33 This also matches the ground state structure factor convergence as G → 0, which has been extensively explored in the QMC and DFT literature. 31,34−36 There is also very strong evidence from the QMC literature that there is a term of N −2/3 in the correlation energy. 30,37,38 The N −2/3 scaling first appears in the exchange energy convergence into the TDL when using an Ewald interaction. 29 To reach an N −1 scaling in the total energy, it is reasonable to assume that the correlation energy must have an equal and opposite term in N −2/3 . For periodic coupled cluster theory, an N −2/3 convergence in the correlation energy would mean that there is an S(G) ∼ G scaling behavior in the transition structure factor that has not yet been identified. This leaves open the question as to whether the CCSD energy has analogous relationships to energies from QMC.
In this study, we will identify how an N −2/3 term arises in the CCSD correlation energy by first analyzing the transition structure factor. We show that the correlation-only transition structure factor (from finite CCSD calculations) fits a functional form with S(G) ∼ G scaling in the limit as G → 0. We then show that the term in G can be canceled by the exchange component of the ground state structure factor, giving rise to the expected overall S(G) ∼ G 2 scaling of the ground state structure factor. We present numerical and analytical results to show how this scaling in S(G) affects energy extrapolations to the thermodynamic limit, paying particular attention to comparing N −1 and N −2/3 extrapolations in practical contexts. We argue in favor of incorporation of an N −2/3 term in the correlation energy extrapolation provided that its prefactor can be derived from exchange-energy calculations.

Coupled Cluster Theory and the Uniform Electron
Gas. All calculations in this paper were performed using coupled cluster theory, where we followed the methods detailed in our previous papers. 19,39,40 Here, we will just outline some of the main methodological details for clarity. In coupled cluster theory, an exponential ansatz is used for the wave function: Ψ = e T̂Φ 0 , where Φ 0 is the ground state wave function, typically taken to be the Hartree−Fock wave function, and T̂is the excitation operator. This wave function is then used to find the coupled cluster correlation energy by projection, i.e., E = ⟨Φ 0 |H|e T̂Φ 0 ⟩. As the work presented here is performed in the uniform electron gas (UEG) where singles excitations are zero due to conservation of momentum, we typically truncate the Tamplitudes to just the doubles to give the coupled cluster doubles (CCD) energy. This energy, then, is calculated using the following equation where t ijab is the T-amplitudes only for the doubles excitations, and v ̅ ijab is the antisymmetrized four-index electron repulsion integrals. Per convention, i and j index occupied orbitals, and a and b index virtual orbitals for a finite basis set. Following a similar derivation to the one in the work by Liao and Gruneis,27 this energy expression is equivalent to one rewritten in terms of the transition structure factor, S(G): The ′ symbol denotes that the sum does not include the G = 0 term. The structure factor is given by . The T appears in place of t to reflect that the indices are now spatial orbitals. The Θ ijab (G) indicates that it only goes over the excitations that are related to the momentum transfers, G, with G being the magnitude of the momentum transfer between the i,j to a,b excitation (i.e., |G| = G). The additional factor of 1/2 comes from the convention we used for the UEG structure factors (for consistency with our prior work). 23 Our electron gas also follows the same setup as described in our previous work, 23,40 with the exception that this work also contains open-shell electron configurations. For our UEG system, we use a simple three-dimensional cubic box with electron numbers that correspond to open-and closed-shell configurations (relative to a grid centered at the Γ-point). The volume of the box is calculated using the Wigner-Seitz radius, r s , such that L r N s , where L is the length of one side of the box. We work exclusively in a plane wave basis set for our UEG calculations, where all the orbitals are described using the relationship k r exp( 1 ) j j · , where k j is a momentum vector for orbital j, and r is the electron coordinate. Ewald interactions are employed for the periodic boundary condition calculations, as per convention, which causes 1/G 2 -type matrix elements to appear in the electron repulsion integrals, v ijab . As there is conservation of momentum in the UEG, only the excitations that correspond to G's that meet the requirement k a − k i = k j − k b = G are used. As with our previous work, we explicitly calculate and include the Madelung term, v M . 23, 29 We also use a finite basis set defined by the M orbitals that lie inside a kinetic energy cutoff E k m cut, . The Hartree−Fock eigenvalues for the occupied and virtual orbitals follow the same conventions as our previous work 23,40 and are lowered in energy by the v M term. In the thermodynamic limit, v M → 0.

Connectivity Twist Averaging.
Twist averaging is typically used to help reduce finite size effects by reducing the fluctuations in the wave function as the system converges to the thermodynamic limit (TDL). 7 Here, the average involves a sum over N s coupled cluster calculations (where N s is the number of twist angles). This increases the cost of running twist-averaged coupled cluster by a factor of N s . In order to help reduce this cost while still obtaining twistaveraged energies for larger systems, we instead use our connectivity twist averaging (cTA) method, which was introduced in other studies. 22,24 With this method, we find a special twist angle for each calculation that reproduces the twistaveraged energies. The method works through evaluating the momentum transfer vectors between the occupied and virtual space, dubbed the "connectivity". These momentum transfer vectors are used to find the twist angle that most closely matches the averaged connectivity using a residual difference calculation. As each system size will have a different connectivity, a special twist angle must be selected individually for each system. The advantage here is that we are now running a calculation using a single twist angle for each system but reproducing twistaveraged energies, lowering the cost of obtaining twist averaged energies by a factor of N s . With this cost reduction, we can then obtain twist-averaged energies for much larger systems, which is vital to our work presented here. 2.
3. An Improved f cut Basis Set Scheme. In our previous work, 19,23 we used a basis set scheme that employs a cutoff factor (f cut ) to truncate the basis set to a given number of orbitals, M. This cutoff factor was chosen such that f cut = E cut,M /E cut,N , where E cut,M refers to the energy cutoff for the basis set with M orbitals, and E cut,N refers to the energy cutoff for the system size containing N electrons. With this method, the E cut,M was calculated manually before being provided to our coupled cluster code for use in truncating the basis set. This basis set scheme will be referred to as f cut (E) in the text, where the (E) is referencing the use of the energy cutoffs to truncate the basis set.
A more precise way to truncate the basis set is to control the number of basis functions per electron with the benefit of allowing the automatic adjustment of the basis set when the electron number changes. In this new basis set scheme, we use our chosen f cut and the number of electrons to calculate M on the fly using the following equation: Here, f cut (M) is the ratio of basis functions per electron rescaled into energy units (using the 3/2 power). With this new method of truncating the basis set, we get a more accurate number of orbitals given our target f cut . One of the advantages of this new f cut (M) basis set scheme is that, unlike with our previous f cut scheme, we are not limited to only the closed-shell system sizes�system sizes determined by symmetry such as N = 14, 38, 54, etc.�as the basis set cutoff is now calculated directly. This allows us to use open-shell systems that break symmetry such as N = 26, 46, 60, etc. Figure 1 shows the results of a comparison between the energies for the f cut basis set. These results support the idea that the f cut (M) basis set helps reduce basis set incompleteness error (BSIE) that causes changes in the finite size error when the electron number changes. Overall, f cut (M) provides a smoother TDL convergence.

Correcting Basis Set Incompleteness Error.
Basis set incompleteness error (BSIE) was handled in the normal way, 23 through deriving a correction to the BSIE from the complete basis set (CBS) limit. This helps ensure that the energies are converged with respect to the basis set before they are extrapolated to the TDL. Since the extrapolations in M and N tend to be independent and commute, 25 we can calculate a basis set correction by choosing a fairly large electron number (here N = 216) and running calculations with increasing basis set sizes. These energies are then extrapolated to the CBS limit, and a basis set correction is calculated using the following equation where M is our chosen basis set size determined by the f cut (M) (here, f cut (M) = 2), and E CBS is the energy at the CBS limit. This correction term, ΔE CBS , is uniformly added to the energies for all N. This process is then repeated for all densities.
2.5. Background Literature on the Ground-State Structure Factor. In our previous work, 23 we fit the CCSD transition structure factor to the following function (inspired by screened MP2) which had a limiting form of S(G) ∼ G 2 as G → 0: This equation leads to an N −1 form for the energy as it approaches the TDL due to the G 2 asymptotic behavior at small G. Here, as we are investigating N −2/3 and N −1 , we require a different functional form that incorporates S(G) ∼ G as G → 0. We take as our inspiration an accurate and well-fitting functional form for the ground-state structure factor, which was proposed by Gori-Giorgi et al. 35 Their functional form incorporated analytical results from the Hartree−Fock approximation (for exchange), the random phase approximation (for the low momentum region), and Quantum Monte Carlo calculations. The ground-state structure factor includes components that correspond to both exchange and correlation structure factors. For our analysis of transition structure factors, we will be using a modified form of this function that only includes the terms coming from Gori-Giorgi et al. correlation structure factor.
Specifically, when we fit our correlation transition structure factor, we used This differs from the structure factor proposed by Gori-Giorgi et al. 35 in that it does not make an attempt to separate the different spin components of the structure factor; instead, we group all of these terms together to simplify the fitting analysis. While those authors also fixed the number of coefficients for each r s value, we found it necessary to include C 6 and C 7 only for r s = 1.0. These were removed for r s = 5.0. Additionally, the inclusion of C 7 was our own addition. We also neglected any treatment of the cusp condition�mainly because the cusp was not our focus in this study, and we used small basis sets for these fits.
In common with Gori-Giorgi et al., 35 the term that is linear in G is constrained such that in the limit of low G, That this is constrained a priori forces the condition that the low-G limit of the exchange-correlation transition structure factor goes as ∼ G 2 . Unlike those authors, we did not fix any further higher-order terms, instead relying upon fitting the function to our data to determine the superlinear coefficients.
We also used an exchange structure factor of the form where q F is the Fermi wave vector equal to q F = α/r s , where α = (9π/4) 1/3 . Here, there are no fit parameters. The momentum transfer G = 2q F is the largest momentum transfer that fits inside the Fermi sphere; thus, we are guaranteed to not have any contributions from G larger than this, allowing S x (G) = 0 for G > 2q F . To make a connection between HF exchange structure factor and transition structure factor, we can note that an exchange transition structure factor can be defined: In this equation, the sum over i and j runs over occupied orbitals, and G is the momentum difference between k i and k j (i.e., the magnitude of G = k i − k j ). The ′ symbol on E x denotes that there is no Madelung contribution to the energy. The ′ symbol on the first sum denotes that the sum does not include i =j and the ′ symbol on the second sum denotes that the sum does not include the G = 0 term. The factor of 1/2 maintains consistency with derivations from other authors and is included due to double counting in sums over electron pairs.

Calculation Details.
For the rest of this work, we will be working with data collected over a range of seven r s values. All r s values were run at an f cut (M) = 2 for a range of N to obtain convergence to the TDL.
The All calculations were performed using a locally modified version of a github repository used in our previous work: http:// github.com/jamesjshepherd/uegccd. 39,40 Hartree atomic units are used throughout.
All graphs were plotted using matplotlib with Python 3.7.3. For the extrapolations to the TDL, the numpy and scipy libraries were used with Python 3.7.3.
All fits for the extrapolation schemes were performed using the curve_fit function from the scipy library in Python. The error in each TDL estimate was calculated from the variance in the fitting parameters.
3.2. Fitting the Transition Structure Factor and Accounting for Exchange. In Figure 2, we show the transition structure factors for r s = 1.0 (N = 508) and another at r s = 5.0 (N = 350). These show calculations of the transition structure factor from CCSD calculations using the relationships described above in eq 1 and eq 2.
In both Figures 2(a) and 2(b), the raw correlation structure factor (S c (G)) data is shown with our transition structure factor fit for both r s = 1.0 and 5.0, respectively. The S c (G) fit is based on a modified form of the ground-state structure factor fit proposed by Gori-Giorgi et al., 35 which only includes the correlation components (see eq 7). Here, the S c (G) fit for both r s has a fixed linear term that is equal and opposite the known linear term from the exchange structure factor. This fixed linear term gives S c (G) a linear convergence to zero as G → 0.0. The close fit between the curve and the data demonstrates that the functional form is consistent with our data. We tested releasing the constraint on the size and sign of the linear term. Unconstrained fits of both result in a curve that is much less well fit, but constraining the sign of the linear term results in a reasonable fit with a coefficient of the same order of magnitude as the original linear term.
Fitting the transition structure factor also allows us to show what happens when the exchange structure factor is included, which is shown in Figures 2(c) and 2(d). Here, the exchange structure factor is plotted using eq 8. These are then combined with the transition structure factor to make the exchangecorrelation transition structure factor: Here, S c (G) is taken from fitting eq 7 to CCSD data, and S x (G) is from the analytical form given in eq 8. When S c (G) and S x (G) are added together, the linear terms cancel by construction (compare eq 7 with eq 8), leaving the exchange-correlation structure factor with a quadratic convergence to zero as G → 0.0. The success of these fits goes a long way in demonstrating that the functions proposed by Gori-Giorgi et al. 35 appropriately model the low G regime of the transition structure factor for the coupled cluster correlation energy.

3.3.
Comparing N −2/3 and N −1 Extrapolations Analytically. With a continuous fit for the structure factor, we are able to examine how finite size effects in the correlation energy converge as the system approaches the TDL. Here, we will follow the same derivation as our previous paper 23 using our analytically derived transition (correlation) structure factor to analyze the FSE for the N −2/3 and N −1 TDL convergence rates. All symbolic manipulations and fits in this section were performed in Mathematica. 48 Using the relationship between the energy and the transition structure factor given in eq 2, we start by integrating over the part of the analytical form of the correlation transition structure factor shown in eq 7, that spans from zero to the minimum G present in our data, G′, to obtain the finite-size error present in the correlation energy: Here, the factor of G 2 comes from the G-space volume element in 3D space. As we are only interested in relative errors, the expression here does not consider any constant prefactors, which are considered in the next section.
It is now possible to estimate the amount of finite size error which is left after extrapolation. Extrapolation consists of fitting the energies to a linear function of the system size (e.g., ). Thus, the removed FSE can be related to the derivative of the function form of the energy (F(G′)) and the energy left after extrapolation is (13) Here, it is important to note that We found analogous expressions for N −1/3 and N −1 . 23 The overall result of this analysis is to be able to test different power laws and how they fit to our structure factor, which was forced to behave as S c ∼ G as G → 0. Figure 3 shows the result of this analysis for the N −1 , N −2/3 , and N −1/3 power laws. The N −2/3 power law shows the best convergence into the TDL, as would be expected given that we fixed the behavior of the structure factor to S c ∼ G as G → 0. While the N −2/3 and N −1 power laws are both reasonably similar, the N −2/3 power law does show slightly improved convergence to the TDL. Here, we can see that the N −1/3 power law overshoots the TDL as the finite size error converges to zero. This is in line with the results from our previous study. 23 This is expected from our use of a structure factor model that goes linearly in G as G → 0.
3.4. Calculating the TDL Energy Using an Interpolation. Given that we now have an analytical form of both the correlation transition structure factor that agree with our current S(G) data as shown in Figure 2, we should be able to integrate over S(G) to get an analytical TDL estimate�in other words, Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article integrating over eq 11 with the constant prefactors included. This is analogous to the approach by Liao and Grueneis. 27 Following this approach, we can calculate the correlation energy as Here, the term is inverse of the k-space volume of one grid point. The additional factor of 1/2 is to maintain consistency with eq 2. The energy produced from this integral, E TDL , is in units of Ha/el. It is important to note here that the E TDL from the above equation still contains basis set incompleteness error, so we need to apply the same uniform basis set correction that we used on the correlation energies. For r s of 1.0 and 5.0, we get TDL energies of −56.52 mHa/el and −22.64 mHa/el, respectively, after all corrections. With these analytical fits, we can now assess which extrapolation scheme from the following section gives the best TDL estimate for our data.

Overview of Extrapolation Schemes to Reach the TDL.
The goal of the rest of this manuscript is to use data from a range of calculations to compare different power laws for their effectiveness at converging exchange and correlation energies to the TDL. From the above analysis, our hypothesis is that N −2/3 is the limiting power law to the TDL, which replaces the N −1 power law that we and other authors have used in recent history for extrapolating the correlation energy. We wish to explore other questions, such as when it is best to use N −2/3 compared with other power laws, and whether extrapolations that have more than one power law are effective.
We will be comparing five different ways to extrapolate the correlation energy to the thermodynamic limit. Here, we give a complete description of how each extrapolation was performed along with a label for each scheme that will be used throughout the rest of this work. We have five schemes in total. Each has a number and may have a letter. The number of the scheme refers to the number of variables used in the fit, while the letter distinguishes different power laws with the same number of variables. For example, Scheme 1A and Scheme 1B both have one variable used in their fit, and they differ because the limiting power law they use is different.
In Scheme 1A, we use a straightforward N −1 convergence rate to extrapolate our basis-set-corrected twist-averaged correlation energies to the TDL. This is the most common way of extrapolating the correlation energy used in the literature, though in some cases has other supporting functions. 4,19,27−31,34,37,41,42,49−60 In all the Scheme 1A extrapolations shown in this work, the extrapolation is performed using the following equation where E TDL (1A) is the energy at the thermodynamic limit for Scheme 1A.
In Scheme 1B, we use an N −2/3 convergence rate to extrapolate to the TDL using our basis-set-corrected twistaveraged energies. This is also a common extrapolation scheme and is often used for the exchange energy. 30,37 For these extrapolations, we used a similar equation to Scheme 1A: The term E TDL (1A) is the energy at the thermodynamic limit for this extrapolation scheme.
In Scheme 2A, we use both the N −1 and N −2/3 convergence rates to extrapolate to the TDL using our basis-set-corrected twist-averaged energies. The equation for the extrapolation then is a combination of Scheme 1A and Scheme 1B: Here, both A and B are free fit parameters that are optimized to give the slopes for the two convergence rates. In theory, the B slope should be similar to the slope of the exchange data, allowing the BN −2/3 term, which ultimately cancels the exchange energy convergence in the total energy. In Scheme 2B, we add a correction term to the N −1 extrapolation from Scheme 1A. The correction term is derived from the twist-averaged exchange energies, which were collected for a range of system size from N = 26 to 826 for r s = 0.1, N = 26 to 508 for r s = 1.0, and N = 26 to 946 for all other r s . The correction term is derived by fitting the exchange to a E x (N) = −B x N −2/3 + E TDL fit. The slope of this fit, B x , is then incorporated into our N −1 extrapolation to the TDL for the correlation energy using the following equation: Here, the B x N −2/3 term is a correction term to help remove some of the residual FSE in the extrapolation and shifts the TDL energy to be more negative. Scheme 3 is an equation with three power laws and is based on the extrapolation scheme presented by Ruggeri et al. 37 In their paper, they suggested the following relationship: We retained the first three terms of this expansion:

Journal of Chemical Theory and Computation
pubs.acs.org/JCTC Article 3.6. Extrapolation Scheme Convergence Across System Size. Having introduced the extrapolation power laws that we would like to test, in this section, we will show the data from our calculations on a variety of UEG systems (Sec. 3.1). The purpose of this section is to compare the extrapolations schemes shown above with one another. In order to compare the effectiveness of extrapolations across different system sizes (modeling an artificial truncation of the data set), we will use a technique we called a windowed extrapolation, where a moving window of points is extrapolated to the TDL using a power law extrapolation scheme.
These windowed extrapolations were performed as follows. Consider N i to be the i th electron number in the data set. If Δi is the number of points in the window to be extrapolated, then the first and second available extrapolations in the data set are over the interval i = [1, 1+Δi] and i = [2, 2+Δi], respectively. The window size, Δi, is typically chosen as the smallest window size that still offers reasonable errors in the fits. The predicted TDL for a single window is assigned to the largest N in the window (i.e., N i+Δi ). For our one-variable extrapolations, Δi = 4 was sufficient. For the two-variable extrapolation, Δi = 6 was used instead. The effectiveness of an extrapolation was then judged by its ability to reproduce the analytically derived energy value calculated in the previous section (Sec. 3.4) and the speed of convergence with system size. We found that these differences in general varied in size between 0.2 mHa/el and 10 mHa/el.
The results of our windowed extrapolations using the extrapolations schemes from Sec. 3.5 are shown in Figure 4 for two densities. Starting with Figure 4(a), where r s = 1.0, we show the convergences of the five extrapolation schemes to the analytical TDL result (from Sec. 3.4). Comparing the different extrapolation schemes, Scheme 1A and Scheme 1B show the slowest convergence. Scheme 1A does not end up agreeing with the TDL within the range of electron numbers we studied, while Scheme 1B only agrees at the largest N. The other three schemes, Scheme 2A, Scheme 2B, and Scheme 3 have a faster convergence, with Scheme 2B and Scheme 3 showing agreement within error to the analytical TDL within the last four points. We also see from this graph a reasonable agreement between the predicted TDL values for Scheme 2B and Scheme 3 across all the windowed extrapolations, indicating that the three terms in Scheme 3 are accounting for the exchange contribution in Scheme 2B. Scheme 2A shows the quickest convergence to the TDL but has oscillatory convergence due to having a free fit on the N −2/3 power law.
We see very similar trends with the r s = 5.0 data shown in Figure 4(b). Here, once again, Scheme 1A and Scheme 1B have the slowest convergence to the analytical TDL compared with Scheme 2B and Scheme 3. In contrast to r s = 1.0 data, however, Scheme 2A shows a closer convergence rate to Scheme 1B here and a wider spread to the TDL values, which results in a slower convergence to the TDL. We also see that the agreement between the Scheme 2B and Scheme 3 extrapolated TDL is maintained with this second r s . This, again, supports the idea that both schemes are accounting for the N −2/3 contribution in the correlation energy.  ) to simulate what would happen if we had less data than we actually have. This allows us to evaluate how each of the extrapolation schemes performs. For each electron number N max , the extrapolation was performed over the Δ i largest system sizes N < N max . The interval Δ i was four for Schemes 1A, 1B, 2B, and 3 and six system sizes for Scheme 2A. Each TDL estimate has been graphed against the largest system size (i.e. N max ) in the range of system sizes used in the extrapolation. This is shown for an r s of (a) 1.   The TDL energies for each electron number N max were obtained using a windowed extrapolation technique (see text). For each N max shown in the table, the difference between TDL energies was taken such that ΔE = E TDL − E Exact , where E TDL is the predicted thermodynamic limit energies at that N max for each extrapolation scheme, and E Exact is the analytical TDL value at each density. For r s = 1.0, the analytical TDL energy is −56.52 mHa/el for the correlation energy and −458.17 mHa/el for the exchange energy. For r s = 5.0, the analytical TDL energy is −22.64 mHa/el for the correlation energy and −91.63 mHa/el for the exchange energy. The analytical TDL energies for the correlation and exchange energies were added together at each r s to get the exchange-correlation analytical TDL value. The number in the parentheses is the error in the difference. All energies are in mHa/el. To make a more detailed comparison, the differences between each of the points for the windowed extrapolations and the analytical TDL values (i.e., ΔE = E TDL − E exact ) are shown in Table 1.
In this table, we show the results of this difference for all extrapolation schemes at r s = 1.0 and r s = 5.0. From these results, we see that Scheme 2B shows the best comparison to the analytical TDL across both r s . Looking at just the r s = 1.0 data, we see that the differences for Scheme 2B shown in the table show a fairly steady convergence of the extrapolated TDL to the analytical TDL as N increases, which is the same as what we saw in Figure 4(a), with agreement within 1 mHa/el (within error) reached by N = 270. Scheme 3 for this r s shows similar trends as Scheme 2B. In contrast, Scheme 1A shows the largest differences to the analytical TDL across all system sizes. For Scheme 1B, the differences show that there is agreement within 1 mHa/el (within error) at system sizes as small as N = 270. Scheme 2A shows 1 mHa/el or less agreement to the analytical TDL at the smallest N and maintains this for all N, but the differences show significant oscillatory behavior as N increases with the largest error seen in the differences across all schemes.
The energy differences from r s = 5.0 show similar results to r s = 1.0, with a few notable differences. Here Scheme 2B and Scheme 3 are slightly different, with Scheme 2B having the smaller differences than Scheme 3 for N ≥ 138. Furthermore, at this density, all of the schemes show differences less than 1 mHa/el by N = 216. Scheme 2B still generally shows the smallest difference out of all five schemes across N starting at N = 174, with the differences dropping to <0.1 mHa/el starting at N = 270. Scheme 2A can also produce energy difference this low but shows nonmonotonic behavior with large errors as N increases that we were seeing at r s = 1.0, making it a less ideal scheme when extrapolating to the TDL. Overall, these results help support the idea that Scheme 2B is the best scheme across densities, with Scheme 3 and Scheme 1B also working well for smaller densities. Interestingly, Scheme 1A, which is the most commonly used extrapolation scheme for the correlation energy, performs the worst out of all the schemes, though at low densities (i.e., r s = 5 or greater) this does not seem to matter as much given that all the extrapolation schemes agree within 1 mHa accuracy.

One-Variable Fits of Exchange, Correlation, and Exchange-Correlation Energies.
There is a tendency for the correlation and exchange energies to mirror one another in how they converge to the TDL. This can be most clearly seen by plotting the two as they converge on the same graph ( Figure 5). We wanted to investigate the hypothesis that, if exchange and correlation energies were extrapolated over the same range and with the same power laws, the error from using Scheme 1A or Scheme 1B would cancel.
Data to investigate this are also shown in Table 1 as exchange and exchange correlation energies. As with the correlation energy, windowed extrapolations were done on both energies with the single-variable power law expansions Scheme 1A and Scheme 1B; the difference to the analytical TDL was taken.
In the case of the exchange energy, we calculated the analytical TDL energy using the known result: 64 E r The exact exchange-correlation energy was then calculated as the correlation energy TDL from eq 14 added to the analytical exchange from eq 21.
Examining the data, we note that the extrapolated exchange TDL underestimates the analytical TDL as convergence is attained. This is in contrast to the fact that the extrapolated correlation TDL systematically overestimates the analytical TDL. These trends were observed for both Scheme 1A and Scheme 1B. The sign difference here corresponds to the mirroring seen in Figure 5.
We also investigated directly extrapolating the exchangecorrelation energy with the one-variable fits, which was added to Table 1. We can see that, overall, this results in Scheme 1A generally giving the better TDL estimate (compared to Scheme 1B) with the smallest residual error. This is consistent with the exchange-correlation energy behaving as N −1 into the TDL. We can also see from this data that there is a similarity between the result from these extrapolations on the exchange-correlation energy and the result of adding together the residual errors after separate extrapolation of the exchange and correlation energies. This suggests that there is a cancellation of error between exchange and correlation energies when they are both extrapolated with the same schemes, such as a N −1 scheme. The convergence of Scheme 1A on the exchange-correlation energy is consistent with the convergence of Scheme 2B on the correlation energy alone.

Correlation Energy TDL Across Densities.
In the previous section, we showed that, out of the five extrapolation schemes, Scheme 2B shows the best comparison to the analytical TDL energies across densities, followed closely by Scheme 3. Here, we want to compare the TDL energies from our extrapolation schemes with the exact TDL correlation energies across various densities. This comparison will give us more evidence for whether or not Scheme 2B is a good generalpurpose extrapolation scheme. For this comparison, we collected TDL predictions across a range of densities (r s = 0.1 to 50.0) using an f cut (M) = 2 for each of our extrapolation schemes described in Sec. 3.5. Calculation details can be found in Sec. 3.1. Exact values come from the Ceperley and Alder results (r s > 1.0) and the Gell-Mann and Brueckner results (r s < 1.0) provided by Perdew and Zunger. 61−63,65 Figure 6 shows the comparison between the different extrapolation schemes and the exact correlation energies across densities.
All schemes are very similar in their TDL extrapolation at r s ≥ 5.0, similar to what we saw in Table 1. This is encouraging, as it means that, for sufficiently high electron numbers, the extrapolations all agree at density ranges that are relevant for everyday materials. For r s < 5.0, the extrapolation methods begin to have different estimates. At these densities, only Scheme 2A, Scheme 2B, and Scheme 3 are shown to be able to capture the TDL energies within the estimated error from extrapolation. All three of these, notably, include a contribution from the power law of N −2/3 . Scheme 2B has a lower error from fitting than Scheme 2A because it has a fixed coefficient derived from exchange energies.
From these results, it would appear as though Scheme 2B is the most consistent in terms of performance when an exchangeenergy-slope can be measured. We emphasize that this is the slope of the exchange energy after a Madelung term has been added. Additionally, in our data set, the exchange energy was also twist-averaged, which may influence the quality of the fits.
3.9. Practical Implications. From our data and analysis, we make the following suggestions: Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article 1. Assuming a situation where more exchange energy data is available than correlation energy data, Scheme 2B is preferred. This is the scheme where the exchange slope is computed separately and then included in the correlation energy extrapolation. Scheme 2B seems especially beneficial at low particle numbers.
2. If there are comparable amounts of both exchange-energy and correlation-energy data, it is advantageous to extrapolate the exchange-correlation energy directly using Scheme 1A. Separately extrapolating exchange and correlation energies using a consistent power law (i.e., using Scheme 1A or Scheme 1B consistently for both exchange and correlation) appears to result in a fortuitous cancellation of error. If extrapolations in the literature were to follow the UEG, therefore, this means current extrapolations are likely to be accurate for the total energy.
3. When the prefactors for both exchange and the leadingorder total energy contributions are known, Scheme 3 can also be used to improve the fit quality significantly over the previous power laws at small electron numbers.

DISCUSSION AND CONCLUDING REMARKS
In conclusion, we incorporated a description of the ground-state (Hartree−Fock) exchange into the transition structure factor of coupled cluster theory. This allowed us to show that there is likely a linear (in G) convergence of the transition structure factor S(G) into the origin, rather than the quadratic (G 2 ) convergence described by previous studies. Using a new basis set cutoff scheme with our previous twist angle selection approach, we calculated unprecedentedly noiseless energy data into the TDL. This allowed us to investigate and compare five schemes for extrapolating the correlation energy into the TDL. We find that some accounting for the N −2/3 term in the extrapolation improves the TDL estimates of the correlation energy. However, we also showed that if the correlation and exchange energies are both consistently extrapolated with an N −1 power law, then the resulting error from using the wrong power law in both cases seems to cancel, at least for the uniform electron gas. As this manuscript was under review, it was also noted to us that an analogous observation of S(G) ∼ G exists in the literature for the random phase approximation, 66 which they formulate as a ring-diagram-based coupled cluster theory. 67 In this work, 66 Bishop and Lurhmann examine the high density electron gas and show analytically that there is a term with appropriate scaling to cancel part of the exchange energy. They further make the identification that the energy density in momentum space (i.e., components of 1/2S(G)v(G) grouped and summed by G) scales linearly in G. Since the number of these terms is proportionate to G 2 and v(G) ∝ 1/G 2 , the result is consistent with our observation of S(G) ∼ G.
We conclude with two limitations of this study. First, as with our previous study, we did not employ any finite size corrections prior to using our extrapolation schemes. We made this choice as we first wanted to see how the extrapolation schemes behaved without adding in additional corrections when accounting for the N −2/3 term in the correlation energy. Second, this study was performed solely on the UEG. It will be important in future studies to show how well these results translate to real materials, including semiconductors and insulators.

■ ASSOCIATED CONTENT Data Availability Statement
The data set used in this work is available at Iowa Research Online at URL: https://doi.org/10.25820/data.006197. The authors declare no competing financial interest. ; B X is found by fitting the exchange energy; and h 2 and t 3 (bolded in the figures for emphasis) are found from external sources. 30,34