Inverse Volume Scaling of Finite-Size Error in Periodic Coupled Cluster Theory

Coupled cluster theory is one of the most popular post-Hartree-Fock methods for ab initio molecular quantum chemistry. The finite-size error of the correlation energy in periodic coupled cluster calculations for three-dimensional insulating systems has been observed to satisfy the inverse volume scaling, even in the absence of any correction schemes. This is surprising, as simpler theories that utilize only a subset of the coupled cluster diagrams exhibit much slower decay of the finite-size error, which scales inversely with the length of the system. In this study, we review the current understanding of finite-size error in quantum chemistry methods for periodic systems. We introduce new tools that elucidate the mechanisms behind this phenomenon in the context of coupled cluster doubles calculations. This reconciles some seemingly paradoxical statements related to finite-size scaling. Our findings also show that singularity subtraction can be a powerful method to effectively reduce finite-size errors in practical quantum chemistry calculations for periodic systems.


I. INTRODUCTION
In the past few decades, ab initio methods for quantum many-body systems, such as density functional theory (DFT), quantum Monte Carlo methods, and quantum chemistry wave function methods, are becoming increasingly accurate and applied to ever larger range of systems [1,2].Unlike molecular systems, periodic systems, including solids and surfaces, require calculating properties in the thermodynamic limit (TDL), a theoretical state in which the system size approaches infinity.However, the TDL cannot be directly accessed in practical applications.Finite-sized computational supercells are employed to approximate this limit, which introduces finite-size errors into the calculations.Finite-size errors can significantly affect the accuracy of calculations, even for systems with thousands of atoms.An extreme case is moiré system such as magic angle twisted bilayer graphene (MATBG), where each computational unit cell consists of approximately 10,000 atoms, and the supercell needs to have more than 100,000 atoms to capture subtle correlation effects [3,4].Directly tackling finite-size effects by enlarging the supercell size is very demanding, even for relatively inexpensive DFT calculations with modern-day supercomputers.For more accurate theories, this task is often computationally intractable.Understanding the finite-size scaling, i.e., the scaling of the finitesize error with respect to the system size, and employing finite-size error correction schemes are, therefore, crucial for obtaining accurate results using moderate-sized calculations.
The sources of finite-size errors in ab initio calculations are multifaceted and complex [5][6][7].These errors are influenced by numerous factors, including system characteristics such as whether it is insulating or metallic, or whether it is a three-dimensional bulk system versus a low-dimensional system.Calculations of electron kinetic energy, electron-ion interaction energy, Hartree energy, Fock exchange energy, and electron correlation energy can all contribute to finite-size errors.The first four types are predominantly single particle in nature, while the electron correlation energy is significantly more complex.To a large extent, electronic correlation is short-ranged, and this characteristic has spurred the development of local correlation methods, whose computational cost may scale linearly with system size.However, for ab initio methods to be accurate, they must also effectively account for van der Waals (vdW) interactions [8].In solids, the cumulative effect of weak van der Waals interactions can become a significant contributor to the energy.The convergence of vdW energy follows an inverse volume scaling, implying that the finite-size error is inversely proportional to the volume of the supercell.The origin of finite-size error is not solely confined integrals [44].This approach provides more accurate estimates and can be more applicable than our previous quadrature analysis based on the Euler-MacLaurin formula for HF and MP2 in [31].In the absence of finitesize correction schemes and assuming exact orbital energies at any k point, the study in Ref. [41] concludes that the finite-size errors of MP3 and CCD both satisfy the inverse length scaling O(N Interestingly, for CCD, earlier numerical calculations did not provide conclusive evidence regarding its finite-size scaling, with different studies suggesting either an inverse volume scaling [15,34] or inverse length scaling [17].More recent calculations demonstrate that the electron correlation energy in periodic coupled cluster calculations should follow an inverse volume scaling, even in the absence of finite-size correction schemes.This observation points to a significant gap in the theoretical understanding of the finitesize error, and prompts the central question we address in this paper: How can we reconcile the following seemingly paradoxical facts?
1. Without finite-size corrections, the finite-size error in CCD exhibits inverse volume scaling.
2. Without finite-size corrections, the finite-size error in MP3 exhibits inverse length scaling.This rate is sharp and cannot be improved.

All MP3 diagrams are encompassed within the CCD formulation.
There are several often-cited physical justifications for expecting that the CCD method, and CC theory more generally, may exhibit superior behavior when applied to periodic systems.One such reason stems from the size extensivity of CC theory.A theory is size extensive when the total energy of two noninteracting identical systems, calculated as a combined system, equates to twice the energy of one system computed independently.Unlike theories such as truncated configuration interaction methods, which are not size extensive, truncated CC theory (such as CCD) possesses this advantageous characteristic.Another reason is that CC theory can be formulated in such a way that it does not explicitly depend on orbital energies.One practical consequence of this is that upon the convergence of the coupled cluster iterations, the Madelung constant correction, often used to reduce finite-size effects in many-body simulations, cancels out naturally.Therefore in this scenario, CC without the Madelung constant correction is equivalent to that with the Madelung constant correction.
We would like to clarify that neither size extensivity nor the cancellation of the Madelung constant alone is sufficient to address the aforementioned question.While size extensivity is indeed a desirable property, many methods such as HF, MP2, and MP3, among others, are all size extensive.In fact, given that periodic systems are infinitely large, size extensivity should be viewed as a necessary condition for the applicability of any numerical method to such systems.However, this property does not provide insight into the convergence rate of the finite-size error.The cancellation of the Madelung constant plays a pivotal role here.Nevertheless, demonstrating that CCD calculation (or even MP3) with the Madelung constant correction exhibits inverse volume scaling in finite-size error is itself a significant challenge.This requires the development of new technical tools not currently available in existing literature.Indeed, the development and application of these tools are the main technical contribution of this paper.
In this paper, we elucidate the origin of the inverse volume scaling behavior.Our analysis consists of two steps.First, we investigate the structure of the CCD amplitude equation.We show that the Madelung constant correction, commonly used to reduce finite-size errors in Fock exchange energy and orbital energies, can also be applied to reduce the finite-size error in ERI contractions within the CCD amplitude equation.We establish a connection between the Madelung constant correction and a quadrature error reduction technique known as the singularity subtraction method [45].By subtracting the leading singular terms from the integrands in the numerical quadratures, the Madelung constant correction reduces the finite-size errors in both the ERI contractions and the orbital energies from O(N − 1 that with the Madelung constant correction, the finite-size errors in CCD(n) and converged CCD calculations satisfy the desired inverse volume scaling.
In the second step of our analysis, we observe that upon convergence of the CCD amplitude equations, the Madelung constant corrections to both orbital energies and ERI contractions perfectly cancel each other out for any finite-sized system.This cancellation ensures that the CCD correlation energy remains the same, regardless of whether the Madelung constant correction is applied.Combining this result with the first step, we conclude that the finite-size error of the correlation energy in converged CCD calculations satisfies the desired inverse volume scaling without the need for any additional correction schemes.However, prior to the convergence of the amplitude equations, this perfect cancellation does not occur, and the finite-size error of CCD(n) calculations remains O(N − 1 3 k ).A similar lack of cancellation occurs when the orbital energies take their exact value at the TDL but the ERI contractions are not corrected, resulting in an O(N − 1 3 k ) finite-size error for both converged CCD and CCD(n) calculations studied in Ref. [41].
To validate our theoretical analysis, we perform CCD calculations on a 3D periodic hydrogen dimer system using the PySCF software package [46].Our numerical results support the conclusions drawn from the theoretical analysis and provide further evidence for the finite-size scaling behavior summarized in Table I.1.
Table I.1.Finite-size scaling of different computational theories with and without corrections to orbital energies and ERI contractions.The first two rows refer to the finite-size scaling of the Hartree-Fock (HF) exchange energy, and the remaining rows refer to the finite-size scaling of the correlation energy (excluding the HF exchange).While we focus on the Madelung constant correction in Ref. [31] and this work, any correction schemes that can reduce the finite-size errors in orbital energies and ERI contractions to O(N −1 k ) can also be applied, and the conclusions drawn here remain valid.In particular, exact values of single particle orbital energies at the TDL satisfy the condition above.N/A means that Madelung constant correction does not apply within the theory.

Theory
Correction to orbital energies

Correction to ERI contractions
The paper is organized as follows.Section 2 introduces background knowledge and basic notations.Section 3 first decomposes the finite-size error in CCD(n) calculations into the errors in four basic components and then describes how the four components contribute to the overall finite-size error with the possible Madelung constant correction.Section 4 explains the key ideas of how the Madelung constant correction can reduce the finite-size error in orbital energies and ERI conttractions from a numerical quadrature perspective.Section 5 illustrates the numerical results that corroborate our error estimate.Lastly, Section 6 discusses the implication of our theoretical study and future directions.

II. BACKGROUND
Consider a unit cell and its Brillouin zone as Ω and Ω * , respectively.Denote the associated real-and reciprocal-space lattices as L and L * .To model such a periodic system, the Brillouin zone Ω * is commonly discretized using a uniform mesh K of size N k (known as the Monkhorst-Pack grid [29]).The orbitals and orbital energies (also called bands and band energies) {ψ nk , ε nk } indexed by orbital indices n and momentum vectors k ∈ K can be solved by the HF method.As a common practice, n ∈ {i, j} refers to an occupied orbital and n ∈ {a, b} refers to an virtual orbital.Throughout this paper, we use the normalized ERI: where q = k 3 −k 1 is the momentum transfer vector, the crystal momentum conservation •r |ψ nk is Fourier representation of the pair product.The primed summation over G means that the possible term with q + G = 0 is excluded in the numerical calculation.Using a finite mesh K, the HF orbital energy without any correction is computed as where Ĥ0 refers to the single-particle component of the many body Hamiltonian.
In this paper, we focus on three-dimensional insulating systems with an indirect gap, i.e., ε aka − ε iki > 0, ∀i, a, k i , k a .To simplify the analysis, we assume that the orbitals are exact at any k point and the number of virtual orbitals are truncated to a finite number.In addition, we assume that the exact orbitals and orbital energies in the TDL are smooth and periodic with respect to their momentum vector index k ∈ Ω * .This assumption is a restriction in our current analysis.For systems free of topological obstructions [48,49], these conditions may be replaced by weaker ones using techniques based on Green's functions, or Hamiltonians defined in the atomic orbital basis instead of the band basis [50].

II.1. CCD theory
Based on the reference HF determinant |Φ , the CCD theory represents the wavefunction as where a † nk and a nk are creation and annihilation operators for ψ nk , T ijab (k i , k j , k a ) (commonly denoted as t aka,bk b iki,jkj in literature) is the normalized CCD double amplitude, and k b ∈ K is uniquely determined using the crystal momentum conservation k ) is defined as the root of a nonlinear amplitude equation that consists of constant, linear and quadratic terms.
In practice, the amplitude equation can be solved using a quasi-Newton method [9,51], which is equivalent to applying fixed point iteration to Here ε N k denotes a diagonal operator with entries , and 1/ε N k gives the diagonal operator with 1/ε N k iki,jkj ,aka,bk b .The operator A N k (T ) is referred to as the ERIcontraction map (see the definition in Appendix A).It consists of contractions between ERIs and T , and does not involve orbital energies.Note that both T and A N k (T ) are tensors indexed by (i, j, a, b) and The CCD correlation energy is then defined as where A TDL can be obtained from A N k by taking 1 The image of A TDL is a set of functions of (k i , k j , k a ) indexed by (i, j, a, b).The CCD correlation energy in the TDL is defined in a similar way as Applying n steps of fixed point iteration over Eq. (3) and Eq. ( 4) with zero initial guess (i.e., t = 0, T = 0), we obtain the CCD(n) amplitude and the CCD(n) energy, CCD(n) is related to the perturbative description of CCD and consists of a subset of finite order perturbation energy terms in the Møller-Plesset perturbation theory.For example, CCD(1) can be identified with MP2, and CCD(2) contains all the terms in MP2 and MP3, and a subset of terms in MP4.
One main result of this paper is the rigorous analysis of the finite-size error in CCD(n) calculation with any fixed n > 0. If the fixed point iterations in both the finite and TDL cases converge to the CCD amplitudes with n → ∞, i.e., T N k n → T N k * and t n → t * (the technical definition of this convergence is provided in Appendix C), the finite-size error analysis for CCD(n) also applies to the converged CCD calculation.In other words, we analyze the finite-size error in CCD calculation using a perturbative approach based on the analysis on CCD(n).

II.2. Madelung constant correction
Ref. [41] shows that the finite-size errors in CCD(n) and CCD both scale as O(N − 1 3 k ) when assuming that exact orbital energies are used in the amplitude equation Eq. ( 3).The same finite-size scaling also appears in Fock exchange energy and occupied orbital energy calculations.One common correction to reduce the finite-size errors in the latter two calculations is to add a Madelung-constant shift [6,47,52] to the Ewald kernel.This shift introduces a correction to all involved ERIs in the calculations as Such a correction is triggered only in ERIs which have fully matched orbital indices, i.e., n 1 = n 3 , n 2 = n 4 , and zero momentum transfer, i.e., k 1 = k 3 .The Madelung constant ξ is defined uniquely by the unit cell and the k-point mesh K as where K q is a uniform mesh that is of the same size as K and contains q = 0 and L K is the real-space lattice associated with the the reciprocal-space lattice q + G with q ∈ K q , G ∈ L. Note that ξ does not vary with respect to σ > 0 [6] and this parameter σ is commonly tuned to control the lattice cutoffs in the summation over L * and L K in Eq. ( 7) when numerically computing ξ.
For finite-size orbital energy calculation in Eq. ( 2), the Madelung constant correction gives For the ERI contractions in A N k (T ), this correction when applied will be triggered in six linear amplitude terms (see Appendix B.2).For example, the 4h2p linear term in Eq. ( 5) can trigger the correction and be modified to Collecting the corrections to all the six terms together, A N k (T ) is modified to In the finite-size CCD(n) and converged CCD calculations, the Madelung constant correction can be applied to the orbital energies, ERI contractions, or both in the amplitude equation Eq. (3).As a result, we may have three correction schemes compared to the standard calculation without any correction.
Particularly, applying the Madelung constant correction to both components gives the amplitude equation It can be easily verified that this amplitude equation has the same solution as Eq. ( 3) without any correction.Its CCD solution is thus identical to the one without correction.However, the associated CCD(n) calculation differs and can be interpreted as solving the original amplitude equation Eq. ( 3) using a quasi-Newton method that has a 2ξ-diagonal shift to the Jacobian matrix.

III. MAIN STATEMENTS
We start our error analysis with the CCD(n) calculation with a fixed number of iterations n > 0 and then generalize the analysis to the fully converged CCD calculation.In CCD(n), the finite-size error is quantified by According to Eq. ( 6), this error can be traced back to two sources: the difference in energy operators, G TDL versus G N k , and the difference in amplitudes, t n versus T N k n .Recall that t n and T N k n are the amplitudes of the system in TDL and finite-size cases, respectively.Let us consider the evaluation map, denoted by M K .This map evaluates a tensor-valued function, initially defined on the product space Ω * × Ω * × Ω * , on a finite-sized grid K × K × K. Consequently, the values of the TDL amplitude t n on this finite-size grid, K × K × K, are given by M K t n and are approximated by the finite-size amplitude T N k n .By applying the triangle inequality, we can decompose the finite-size error into two distinct sources: the errors arising from the discretized energy calculation using the exact amplitude M K t n , and the errors stemming from the amplitude calculation itself.
Here we use the fact that |W ijab (k i , k j , k a )| can be upper bounded uniformly by a constant.
To further break down the error in amplitude calculation, we note that t n and T N k n are recursively constructed by Eq. ( 6) with initial values t 0 = 0 and T N k 0 = 0.As a result, the error in the CCD(n) amplitude calculation can also be recursively decomposed using the same strategy above as The three error terms from this dissection can be interpreted as the errors in ERI contractions using exact CCD(n − 1) amplitudes, orbital energies, and CCD(n − 1) amplitude calculation composed with A N k , respectively.For the last term, it can be shown that entries in have the same scaling with respect to N k as those in applying the dissection recursively, we find that the error in the CCD(n) amplitude calculation is determined by those in the ERI contractions and orbital energies, i.e., the first two terms in Eq. (13).
Overall, the finite-size error of CCD(n) calculation can be decomposed into errors in three basic factors: (1) energy calculation using exact CCD(n) amplitude, (2) ERI contractions using exact CCD(n − 1) amplitude, and (3) orbital energies.This error decomposition is also valid when applying the Madelung constant correction to orbital energies (Eq.( 8)) or ERI contractions (Eq.( 10)).By analyzing these three error sources with or without corrections separately, we can obtain the finite-size error estimate for CCD(n) with various correction schemes.
The Madelung constant correction can reduce the finite-size error in orbital energies from O(N This correction is at the HF level.One main technical result of this work is to show that the Madelung constant correction can also reduce the finite-size error in most (but not all) entries of the ERI contraction (note that Now what happens to the finite-size error of the CCD calculation as n → ∞?For gapless and small-gap systems, it has been observed in practice that the fixed point iteration might not converge or the amplitude equation might have multiple solutions.While we focus on systems with an indirect gap, it is worth noting that even in such favorable scenarios, existence and uniqueness of solutions to the CCD amplitude equations in both finite and TDL cases remain an open question and are beyond the scope of this paper. To study the finite-size error of CCD via the above results on CCD(n), we make additional technical assumptions (see Appendix C) that can guarantee the convergence of CCD(n) to CCD.Under these assumptions, we show that the finite-size scaling of CCD calculation is upper bounded by those of its associated converging CCD(n) calculations.Numerical observations (see Section V) further show that this finite-size scaling estimate through CCD(n) is asymptotically sharp for CCD calculation with the Madelung constant correction applied to orbital energies, ERI contractions, or both.

each of the following scenarios (1) the Madelung constant correction is only applied to the ERI contraction A N k (2) the Madelung constant correction is only applied to the orbital energy
nk in the CCD calculation, the overall finite-size error scales as O(N −1 k ).Compared to Theorem 1, Theorem 2 does not address the finite-size error of the CCD calculation when no finite-size correction is applied.The natural conclusion from Theorem 1 is that this finite-size scaling should be O(N − 1 CCD calculation without any correction is equivalent to the CCD calculation with the Madelung constant correction applied to both ERI contractions and orbital energies.Specifically, when applying the Madelung constant corrections Eq. ( 8) and Eq. ( 10), the CCD amplitude equation in Eq. ( 11) can be formulated as With 2ξT on both sides cancelling each other, this reformulation is exactly reduced to the original amplitude equation Eq. ( 3) without corrections.Therefore, the roots of the two amplitude equations with and without the corrections are the same, and the correlation energy in converged CCD calculation without any corrections is the same as that with the Madelung constant correction.In other words, when investigating the finite-size error of CCD calculation without corrections, we should apply Theorem 2 with the Madelung constant correction applied to both A N k and ε N k nk .This yields a sharp estimate O(N −1 k ) and explains the origin of the inverse volume scaling of the finite-size error.
Corollary 3.Under the same additional conditions as in Theorem 2, the finite-size error of the CCD correlation energy without finite-size correction scales as O(N −1 k ).We provide the proof of Theorem 2 in Appendix C, and Corollary 3 follows directly from Theorem 2 and the reasonings above.

IV. KEY IDEAS
As demonstrated in the previous section, the finite-size errors in CCD(n) and CCD calculations can be reduced to the errors in three simpler basic calculations: energy calculation using exact amplitudes, ERI contraction using exact amplitudes, and orbital energies.A key observation is that the finite-size errors in all the three calculations can be interpreted and analyzed from a numerical quadrature perspective.
For a function g over a hypercube V , we denote a (generalized) trapezoidal rule using a uniform mesh X of V and its quadrature error as Under the assumption of exact orbitals at any k point in finite-size calculation, the errors in (1) orbital energy, (2) energy, and (3) ERI contraction A N k can be respectively formulated by their definitions as For ERI contraction in the last line, the error is denoted by a tensor with indices (i, j, a, b, k i , k j , k a ), and the error for the 4h2p linear term Eq. ( 5) is detailed with k k being the integration variable.Meanwhile, terms not shown account for errors from computing other linear and quadratic amplitude terms.From a numerical quadrature perspective, all the three finite-sized calculations approximate the corresponding integrals in the TDL using the trapezoidal rule and a finite mesh K to discretize Ω * .Therefore, their finite-size errors can be estimated systematically by quadrature error analysis.
In general, the quadrature error associated with a trapezoidal rule is influenced by the integrand's smoothness and boundary conditions.If we take ∆h as the mesh size along each dimension, the quadrature error for a smooth integrand is typically of order O(∆h 2 ).Interestingly, if the integrand is also periodic, the error diminishes much more rapidly than it does in a nonperiodic scenario.The decay rate is faster than any finite power of ∆h, showcasing a super-algebraic decay [53].However, for an integrand periodic but marked by singularities, its quadrature error tends to taper off at a slower rate.
For quadrature errors in Eq. ( 14), the involved integrands are all periodic across their integration domains.However, many of these integrands have singularities within the domains, affecting the scaling of their quadrature errors.As to be demonstrated next, these integrands include point singularities arising from both ERIs and amplitudes, resulting in low-order power-law decay of the corresponding quadrature errors as N k increases.

IV.1. Singularity structure and quadrature error estimate
All the quadrature errors in Eq. ( 14) have integrands comprising of either ERIs or contractions between ERIs and exact CCD(n) amplitudes.The integration variables are momentum vectors sampled in the Brillouin zone Ω * .Consequently, understanding the singularity structure of ERIs and exact CCD(n) amplitudes with respect to their momentum vector indices is crucial in comprehending the singularity structure of these integrands and ultimately estimating the quadrature error in the three basic calculations.
First consider a generic ERI ) and treat it as a function of k 1 , k 2 , and q = k 3 − k 1 in Ω * .By its definition, the ERI can be separated as Since we assume all orbitals ψ nk periodic and smooth with respect to k ∈ Ω * , this ERI has a point singularity at q = 0 only due to the first fraction term.Specifically, this term is of fractional form f (k 1 , k 2 , q)/|q| 2 with a smooth numerator f .As can be verified by direct calculation, such a fraction term can have its point singularity at q = 0 described by the following general concept called algebraic singularity.
where α denotes a non-negative d-dimensional derivative multi-index.For brevity, f is said to be singular at x 0 with order γ.See Definition 8 in Appendix for the rigorous mathematical definition.
A representative example of such a singular function with order γ is g(x)/|x| 2 where g(x) is smooth and scale as O(|x| 2+γ ) near x = 0. Using orbital orthogonality, the generic ERI exhibits singularities at q = 0 with order 0 when band indices mismatch (n We can now characterize the singularity structure of integrand i W inin (k i , k, k i ) defined by orbital energy ε TDL nk in Eq. ( 14).Fixing (n, k), the leading singularity in each W inin (k i , k, k i ) (as a function of k i ) comes from the exchange term ik i , nk|nk, ik i which has algebraic singularity at k i = k with order 0 when i = n and −2 when i = n.In computing ε TDL nk , the overall integrand is thus singular at k i = k with order 0 for a virtual band n and −2 for an occupied band n.
For such periodic functions with one point of algebraic singularity, the conventional textbook analysis of the trapezoidal rule is overly pessimistic.A key technical tool in this work is Lemma 5 below, which provides a rigorous and sharp quadrature error estimate linking its error scaling to the singularity order.(See Lemma 22 in Appendix for a more general statement.) ] d and smooth everywhere except at x = 0 with order γ −d + 1.At x = 0, f (x) is set to 0. The quadrature error of a trapezoidal rule using an m d -sized uniform mesh X that contains x = 0 can be estimated as Combining this error estimate with the integrand singularity structure for orbital energies (with m = N 1 3 k , d = 3, and γ = −2, 0 for occupied and virtual orbitals, respectively), we obtain the finite-size error estimate for orbital energy calculation as To adapt the above approach for examining energy and ERI contraction calculations with exact amplitude, we need to characterize the singularity structure of the exact CCD(n) amplitude entries.First we note that the exact CCD(1) amplitude is just the MP2 amplitude As a function of (k i , k j , k a ), each CCD(1) amplitude entry indexed by (i, j, a, b) has the same singularity structure as the included ERI term which is singular at k a − k i = 0 with order 0. It turns out that the exact CCD(n) amplitude with any n > 0 all has the same singularity structure as ak a , bk b |ik i , jk j .
Lemma 6 (Singularity structure of the amplitude, Lemma 4 in [41]).In CCD(n) calculation with n > 0, each entry of the exact double amplitude t n = {[t n ] ijab (k i , k j , k a )} belongs to the following function space f is smooth everywhere except at k a = k i with algebraic singularity of order 0, Combining the above singularity structure characterizations of ERIs and exact CCD(n) amplitudes, we are able to analyze the integrands for energy and ERI contraction calculations in Eq. ( 14).Take the CCD(n) exchange energy term as an example whose finite-size error can be formulated as with integration variables (k i , k j , k a ).For each set of (i, j, a, b), both the associated ERI and exact amplitude in the integrand exhibit an algebraic singularity of order 0 at k a −k j = 0 and k a −k i = 0, respectively.For such a product of two functions with algebraic singularities, we also provide a rigorous quadrature error estimate similar to Lemma 5. Specifically for all the integrands in the energy and ERI contraction calculations, their quadrature errors are determined by the most singular product components in the integrand.The quadrature error scaling still scales as m −(d+γ) similar to Lemma 5 but with γ denoting the minimum algebraic singularity order of all ERIs and exact amplitudes.
For the exchange term above, all involved ERIs and exact amplitudes have point singularities of order 0, and similar for the direct term.Thus, we can get the finite-size error estimate for the energy calculation using exact CCD(n) amplitude as Similar analysis is also applicable to ERI contractions using exact amplitudes.The key distinction lies in the fact that the ERI contractions involve integrands formed by ERIs with stronger singularities, such as those defined by particle-particle or hole-hole diagrams.A prominent example is the 4h2p linear term in Eq. ( 5), where the involved ERI kk k , lk l |ik i , jk j exhibits an algebraic singularity of order −2 at k k = k i when k = i and l = j.Consequently, as per the above analysis, the finite-size error in the 4h2p linear term calculation alone scales as O(N ).This error turns out to dominate the overall finite-size error in the ERI contraction calculation, and we have

IV.2. Madelung constant correction as a quadrature error reduction method
To reduce the quadrature error for a singular integrand, one common numerical quadrature technique is the singularity subtraction method [45].Essentially, this method involves constructing an auxiliary function h that possesses the same leading singularity as the integrand g.The integral is then approximated as: This approximation consists of the numerical quadrature of g − h and the exact integral of h (which can be computed analytically or numerically with high precision).It is also equivalent to adding a correction E V (h, X ) to the numerical quadrature of g.By this correction, the quadrature error changes from E V (g, X ) to E V (g − h, X ).Since h removes the leading singularity of g in the subtraction g − h, E V (g − h, X ) can be asymptotically smaller than E V (g, X ).
The Madelung constant defined in Eq. ( 7) can be reformulated using an arbitrary σ > 0 as Compared to Eq. ( 19), this representation connects ξ to the singularity subtraction correction defined by an auxiliary function The effectiveness of the Madelung constant correction to reduce the finite-size error can be rigorously explained by this connection.
Taking the occupied orbital energy ε nk as an example, its exchange portion leads to the dominant finitesize error and the associated Madelung constant correction modifies the calculation as (with a change of variable Comparing this calculation with Eq. ( 19) and Eq. ( 20), the Madelung constant correction exactly uses the auxiliary function h σ (q) to remove the leading singularity (i.e., 4π/|Ω||q| −2 ) of the target integrand, and thus reduces the associated finite-size error asymptotically to O(N −1 k ).One major technical contribution of this paper is to rigorously prove the effectiveness of the Madelung constant correction for reducing the finite-size error in the ERI contractions calculations following the same singularity subtraction interpretation.For instance, consider the 4h2p linear term calculation in Eq. ( 5) with any fixed entry index (i, j, a, b, k i , k j , k a ).Using the change of variable k k → k i − q, this term with the Madelung constant correction can be detailed as The leading singularity of the integrand comes from the product with (k, l) = (i, j).In this product, the ERI is singular at q = 0 with order −2 and the amplitude is singular at q = k i − k a with order 0. The correction in Eq. ( 21) defines a singularity subtraction with auxiliary function h σ (q)t ijab (k i , k j , k a ).This auxiliary function shares exactly the same leading singularity as the integrand at q = 0 due to the ERI term, i.e., 4π |Ω| Similar to the orbital energy analysis, the finite-size error in this 4h2p linear entry can be reduced to O(N −1 k ).However, the key difference here is that this error reduction is the case only for most but not all the amplitude entries.The exception is when the amplitude singularity q = k a − k i is close or equal to the ERI singularity q = 0, i.e., when computing an ERI contraction entry whose momentum vector indices k i and k a are close or identical.In the worst case when k i = k a , the two singularities overlap and the finite-size error of such 4h2p linear entry can be shown to be still of scale O(N ).Similar analysis also applies to other terms in the ERI contraction calculation.
To summarize, the Madelung constant correction does not uniformly reduce the finite-size errors in the ERI contraction tensor.This is the case for the 4h2p linear term and also for other terms in the ERI contraction.More precisely, we have the following technical error estimate (see Lemma 11 in Appendix for the general statement and proof).
Lemma 7 (Error in ERI contractions).The finite-size error in the ERI contractions using exact CCD(n) amplitude t n with the Madelung constant correction in Eq. ( 10) satisfies As a result of this nonuniform error reduction, the maximum entrywise finite-size error in the CCD(n) k ).However, the average entrywise error satisfies the bound From the error decomposition in Eq. ( 12), such a refined bound is sufficient for our finite-size error analysis.Therefore the Madelung constant correction to the ERI contraction and the orbital energies can effectively reduce the finite-size error in the overall CCD(n) energy calculation to O(N −1 k ).

V. NUMERICAL EXAMPLES
To validate the above theoretical analysis, we conduct CCD and CCD(n) calculations on a 3D periodic system of hydrogen dimers.One hydrogen dimer is positioned at the center of each cubic unit cell with an edge length of 6 Bohr in the x-direction, and a separation distance of 1.8 Bohr.For each uniform mesh K, we perform an HF calculation on K to obtain the orbitals and orbital energies.We then perform CCD and CCD(n) calculations under four distinct settings, each with a different combination of the Madelung constant correction to the orbital energies and the ERI contractions in the amplitude equation.Specifically, we compute CCD and CCD(n) with both corrections, the correction to the orbital energies alone,the correction to the ERI contractions alone, and without any corrections.All the calculations are carried out using the PySCF package [46] with a minimal basis set gth-szv.
Figure V.1 illustrates the numerical results of the CCD(1), CCD(2), CCD(3), and converged CCD calculations.For CCD (1), which is identical to MP2, we have T 0 = 0 and the Madelung constant correction to the ERI contractions has no effect.As a result, the curves for the calculations with and without this correction are identical.For CCD(2) and higher, the four correction settings produce distinct curves.Only the CCD(n) calculation with both corrections exhibits convergence rate of O(N −1 k ), while the other three calculations have convergence rates of O(N − 1 3 k ).This highlights the importance of taking into account the Madelung constant correction to both the orbital energies and ERI contractions in higher-level CCD(n) calculations.As CCD(n) converges to CCD, the difference between the calculation with both corrections and the one without any corrections gradually diminishes to zero, and the finite-size error satisfies the inverse volume scaling.

VI. DISCUSSION
Recent years have witnessed significant progresses in leveraging wavefunction methods to study solids.One important driver of this trend is the potential of these methods to provide systematically improvable results.Among the simplest post-HF methods, MP2 has been applied to increasingly large periodic systems along with the development of efficient implementations utilizing density fitting and parallelization techniques [16,23].Notably, the next order method, MP3, has also been applied to periodic systems recently for the first time [24].Though computationally much more expensive, CC theory has also gained momentum in this context.It is being increasingly employed to compute ground state and band structures for a variety of solid materials including insulators and metallic systems [17-19, 22, 25].Furthermore, related theories such as equation-of-motion coupled-cluster theory, GW method, and algebraic diagrammatic construction theory also start their utilizations in computing excited state properties of solids [20,21,24].Due to the omnipresence of finite-size effects, as wavefunction methods continue to be increasingly applied to periodic systems, the importance of developing a thorough understanding of finite-size effects becomes even more pronounced.
In this work, we fill a gap in understanding the finite-size error in periodic coupled cluster calculation for insulating systems.Notably, we reveal an unexpected inverse volume scaling of this error in CCD calculation.This behavior manifests even in the absence of any finite-size correction schemes, owing to an error cancellation.Our findings, together with the methodologies employed in this study, provide valuable insights for practitioners, method developers, and theorists.
For practitioners, when applying computational quantum chemistry methods to periodic systems, reducing finite-size errors using techniques such as power-law extrapolation requires an in-depth understanding of the error scaling.This is particularly important when calculations are constrained to small-sized systems due to the steep increase of the computational cost with respect to the system size and limited resources.Many production-level quantum chemistry packages that support periodic systems use certain finite-size corrections to reduce errors in the HF exchange energy calculation and the HF orbital energies.For instance, the truncated Coulomb correction scheme [54,55] can be applied to insulating systems, so that the finite-size error in orbital energies decays super-algebraically with respect to N k .However, our analysis shows that if this correction is only applied to the orbital energy but without any correction to ERI contractions, the finite-size error of the correlation energy deteriorates from inverse volume scaling to inverse length scaling.This is also a finding consistent with numerical observations.
For method developers, a key aspect of our result is the connection between the Madelung constant correction and the singularity subtraction method.This relationship serves not just as a crucial element in our theoretical proof, but also points towards new methods for further finite-size error reduction.The Madelung constant correction operates as a one-shot correction, while the singularity subtraction method can be systematically improved to reduce the finite-size error.Our analysis and correction schemes can also be extended to more advanced coupled cluster theories such as coupled cluster singles and doubles (CCSD), and coupled cluster singles, doubles and triples (CCSDT).Exploring the finite-size error in lower dimensional systems, such as magic angle twisted bilayer graphene (MATBG) that requires a different form of the Coulomb kernel [4,56], can lead to different finite-size scaling patterns and novel correction schemes.Going beyond finite-size corrections, the singularity structure of the amplitude, as outlined Lemma 6, may be of independent interest.The singularity structure imposes important analytic constraints that should be taken into account when developing efficient numerical methods for compressing the CC amplitude tensor.
For theorists, several critical questions persist: How can finite-size error analysis be integrated with the study of complete basis limits to tackle basis set dependence? How should the finite-size error behavior in metallic systems be analyzed, especially when the orbital energy difference in the denominator could vanish?Can the scope of finite-size analysis be broadened to include more complex systems like disordered systems and finite-temperature alloys?These questions present a fertile ground for future research.

Appendix A: CCD amplitude equation
In this section, we use a capital letter to denote an index pair consisting of the orbital index and the kindex.For instance, I = (i, k i ), J = (j, k j ), A = (a, k a ) etc. We use P ∈ {I, J, K, L} to refer to occupied orbitals (a.k.a., holes) and P ∈ {A, B, C, D} to refer to unoccupied orbitals (a.k.a., particles).Any summation P refers to summing over all occupied or virtual orbital indices p and all momentum vectors k p ∈ K while the crystal momentum conservation is enforced according to the summand.This notation is used only in this section to simplify the notation and also to connect the equations to those in the molecular case [9] for better readability.
Using a finite mesh K of size N k , the normalized CCD amplitude is defined as the solution of the amplitude equation where , and P is a permutation operator defined as This reformulation of CCD amplitude equation is derived from the CCSD amplitude equation in [30] by removing all the terms related to single amplitudes and normalizing the involved ERIs and amplitudes (which gives the extra 1/N k factor in the equation and the intermediate blocks).The intermediate blocks in the equation are defined as and their momentum vector indices also assume the crystal momentum conservation In the TDL, the amplitude equation for the exact amplitude t * = {t ijab (k i , k j , k a )} := {t AB IJ } as functions of k i , k j , k a ∈ Ω * can be formulated by letting K in Eq. (A1) converge to Ω * as where the intermediate blocks in the TDL are defined as Appendix B: Proof of Theorem 1 Restatement of Theorem 1.In CCD(n) calculation, the finite-size error in the correlation energy scales as O(N k ) in each of the following scenarios (1) there is no finite-size correction, (2) the Madelung constant correction is only applied to the ERI contraction A N k , and (3) the Madelung constant correction is only applied to the orbital energy ε N k nk .When the Madelung constant correction is applied to both A N k and ε N k nk in the CCD(n) calculation, the overall finite-size error scales as O(N −1 k ).As a special case, the same conclusion applies to MP3 calculations.

B.1. Proof Outline
The main context of this paper has provided a brief description of the main proof idea.In this proof, we will recap some of the equations and concepts discussed before to make it self-contained.The proofs of Theorem 1 for CCD(n) calculations with various types of corrections are based on the error splitting in Eq. ( 12) and Eq. ( 13), i.e., Note that the amplitude error here is measured in the average norm as The finite-size error in CCD(n) calculation is thus decomposed into the error in energy calculation using exact amplitude, the error in ERI contractions, the error in orbital energies, and the error accumulated from previous iteration.The latter three errors together make up of the error in amplitude calculation.The error in energy calculation using exact amplitude is studied previously in Ref. [41].For completeness, we provide a brief review of the main results below.

Brief review of error in energy calculation with exact amplitude
One basic observation is that this error in CCD(n) calculation can be interpreted as the quadrature error of a specific trapezoidal rule as Since both W ijab and [t n ] ijab are periodic with respect to k i , k j , k a ∈ Ω * , the asymptotic scaling of this quadarture error depends on the smoothness of these two components that constitute the integrand.A general ERI n 1 k 1 , n 2 k 2 |n 3 k 3 , n 4 k 4 can be viewed as a function of momentum vectors k 1 , k 2 and its momentum transfer vector q = k 3 − k 1 .This function is periodic with respect to each variable in Ω * and its definition in Eq. ( 1) can be decomposed as The numerators of all the fractions are smooth with respect to k 1 , k 2 , q (note the assumption that ψ nk (r) is smooth with respect to k).Therefore, this ERI example is smooth with respect to k 1 , k 2 and has one point singularity in Ω * with respect to q at q = 0, which is due to the first fraction term.For any fixed k i , k j , this point singularity can be characterized using the concept of algebraic singularity of certain orders.
Definition 8 (Algebraic singularity for univariate functions).A function f (x) has algebraic singularity of order γ ∈ R at x 0 ∈ R d if there exists δ > 0 such that where the constant C α,δ depends on δ and the non-negative d-dimensional derivative multi-index α.For brevity, f is also said to be singular at x 0 of order γ.
The numerator of the first fraction in Eq. (B4) scales as O(|q| s ) with s ∈ {0, 1, 2} near q = 0 using the orbital orthogonality.The value of s depends on the relation between orbital indices (n 1 , n 2 ) and (n 3 , n 4 ).As a result, the algebraic singularity of the ERI example above at q = 0 with any fixed k 1 , k 2 has order in {−2, −1, 0}.In addition, to connect this singularity with varying k 1 , k 2 ∈ Ω * , we introduce the algebraic singularity with respect to one variable for a multivariate function.
Definition 9 (Algebraic singularity for multivariate functions).A function f (x, y) is smooth with respect to y ∈ V Y ⊂ R dy for any fixed x and has algebraic singularity of order γ with respect to x at x 0 ∈ R dx if there exists δ > 0 such that where constant C α,β,δ depends on δ, α and β.Compared to the univariate case in Definition 4, the key additions are the shared algebraic singularity of partial derivatives over y at x = x 0 of order γ and the independence of By this definition, the ERI example is smooth everywhere with respect to k i , k j , q ∈ Ω * except at q = 0 of order γ ∈ {−2, −1, 0}.Specifically, the order γ equals to −2, −1 and 0, respectively, when the orbital indices are fully matched, i.e., n 1 = n 3 , n 2 = n 4 , partially matched, i.e., n 1 = n 3 , n 2 = n 4 or n 1 = n 3 , n 2 = n 4 , and not matched, i.e., n 1 = n 3 , n 2 = n 4 .If treating the ERI example as a function of k 1 , k 2 , k 3 instead, we equivalently claim that the function is singular at k 1 = k 3 of order γ.
One key result in [41] is the singularity structure characterization for the exact CCD(n) amplitude t n , which is essential for estimating the quadrature error in Eq. (B3).It turns out that that each exact amplitude entry [t n ] ijab (k i , k j , k a ) with any n > 0 has one point of algebraic singularity of order 0 at k a − k i = 0, sharing a similar singularity structure as the ERI ik i , jk j |ak a , bk b or the exact MP2/CCD(1) amplitude entry (ε TDL iki,jkj ,aka,bk b ) −1 ak a , bk b |ik i , jk j .
Restatement of Lemma 6 (Singularity structure of the amplitude, Lemma 4 in [41]).In CCD(n) calculation with n > 0, each entry of the exact double amplitude t n belongs to the following function space f is smooth everywhere except at k a = k i with algebraic singularity of order 0, f is smooth with respect to k i , k j at the singularity k a = k i .
Based on the above singularity structures of ERIs and exact amplitudes, the integrand in the energy calculation in Eq. (B3) consists of products of periodic functions where each has one point singularity of order 0. (Recall that W ijab is the antisymmetrized ERI and consists of two ERIs that can be treated separately.)We provide a sharp quadrature error bound for trapezoidal rules over periodic functions in such a product form, and its application to Eq. (B3) gives Lemma 10 (Energy error with exact amplitude, Lemma 5 in [41]).In CCD(n) calculation with any n > 0, the finite-size error in the energy calculation using exact amplitude t n can be estimated as where constant C depends on t n .

Error in amplitude calculation
Based on the error splitting in Eq. (B2), analysis of the error in amplitude calculation is reduced to estimating the two main error terms

and understanding how
n−1 from the previous iteration.Like the energy calculation, the errors in ERI contractions and orbital energies consist of specific quadrature errors.In both cases, we can show that the Madelung constant correction is connected to certain singularity subtraction methods and can significantly reduce the corresponding dominant quadrature errors.
Lemma 11 (Error in ERI contractions).In CCD(n) calculation with any n > 0, the finite-size error in the ERI contractions using exact amplitude t n without any corrections has its entries bounded as The Madelung constant correction reduces this error as where q = k a − k i + G 0 with G 0 ∈ L * chosen such that q ∈ Ω * .In both cases, constant C depends on t n but not on the entry index (i, j, a, b) and Remark 12.The prefactor 1/|q| 2 in the above estimate is important when q ∈ K q is near the origin.For

Lemma 13 (Error in orbital energies). The finite-size error in orbital energies without any corrections is bounded as
The Madelung constant correction reduces this error as In both cases, constant C is independent of the entry index n, k ∈ K.
Note that there are three multipliers for the three error terms in the amplitude error splitting Eq. (B2).These prefactors are bounded by constants independent of N k .First, the orbital energy difference in Eq. (B2) satisfies |ε TDL iki,jkj ,aka,bk b | 2ε g by the assumed indirect gap ε TDL aka − ε TDL iki ε g > 0. Second, since A N k (T ) consists of constant, linear and quadratic terms of T , a straightforward estimate shows (e.g., using Lemma 7 in [41]) that max ijab,ki ,kj,ka∈K where the L ∞ -norm of [t n ] ijab is finite according to Lemma 6.
Based on the above estimates of these multipliers, we have that the summation of the first two error terms in Eq. (B2) is dominated by the error in the ERI contractions and the orbital energies discussed in Lemma 11 and Lemma 13.If applying the Madelung constant correction to both orbital energies and ERI contractions, the summation of the first two error terms has its entries bounded asymptotically the same as Eq.(B6).Otherwise, its entries are bounded asymptotically as in Eq. (B5).Lastly, for the third error term in Eq. ( 13), the application of A N k /A N k ,ξ can be proved to maintain the entrywise error scaling obtained in Lemma 11.
Lemma 14.Consider two arbitrary bounded amplitude tensors T, S ∈ C nocc×nocc×nvir×nvir×N k ×N k ×N k and assume their entrywise upperbound independent of N k .If the difference between T and S satisfies an entry bound like Eq. (B5), i.e.,

The ERI-contraction map without any corrections
If the difference between T and S satisfies an estimate similar to that in Eq. (B6), i.e., The ERI-contraction map with Madelung constant correction A N k ,ξ satisfies Combining the estimates of the three error terms in Eq. ( 13) and the initial condition M K t 0 = T N k 0 = 0 in both cases with and without Madelung constant correction, we can obtain the entrywise estimate of the error in the CCD(n) amplitude calculation recursively.First, for the CCD(n) calculation without any corrections or with partial Madelung constant corrections to either orbital energies or ERI contractions, we have Accordingly, the average entrywise error can be estimated as Secondly, for the CCD(n) calculation with the Madelung constant correction to both orbital energies and ERI contractions, we have Accordingly, the average entrywise error can be estimated as Plugging the above estimate and the error estimate for energy calculation in Lemma 10 into Eq.( 12) then finishes the proof of Theorem 1.
Remark 15.In our previous work [41], we use the maximum entrywise norm to characterize the finite-size error in the amplitude calculation, which loosens the average entrywise norm used in the error splitting in Eq. (B1) as Without the Madelung constant correction, the maximum norm suffices since all entries in the amplitude error are of the same scale, as shown in Eq. (B7).However, this norm is no longer sufficient for the calculation with the Madelung constant correction, bounded in Eq. (B8).Because the maximum entrywise error is now of scale O(N k ), while most entries are of scale O(N −1 k ).In this case, the average entrywise norm provides a necessary and tighter estimate of amplitude error.

B.2. Proof of Lemma 11: error in ERI contractions
According to the singularity structure of exact CCD(n) amplitude in Lemma 6, we consider the ERIcontraction using an arbitrary exact amplitude t ∈ T(Ω * ) nocc×nocc×nvir×nvir .Fixing a set of entry index (i, j, a, b) and (k i , k j , k a ) ∈ K × K × K, the error in the indexed ERI-contraction entry can be detailed as (by comparing Eq. (A1) and Eq.(A2)) where the constant terms cancel with each other, the first term is the error in the 4h2p linear term calculation with the Madelung constant correction, and the second term is the error in the 4h2p quadratic term calculation.The neglected ones are the errors in remaining linear and quadratic terms calculations, which can all be similarly formulated as quadrature errors of specific trapezoidal rules.
In the error analysis for CCD calculation without corrections in ERI contractions [41], it has been shown that without the Madelung constant correction the error entry in the ERI contractions is uniformly bounded as More specifically, all the quadratic terms and part of the linear terms that only contain ERIs with mismatched orbital indices contribute at most O(N −1 k ) errors in both Eq.(B9) and Eq.(B10).The dominant error in Eq. (B10) comes from the calculation of the following six linear terms which contain ERIs with fully or partially matched orbital indices.The Madelung constant correction in Eq. ( B9) is exactly triggered in these six terms.In this proof, we focus on the error estimate for the 4h2p linear term (the first term above) with the correction, and similar analysis can be applied to all the other five terms.
Denote the ERI-amplitude product in the 4h2p linear term with orbital indices (k, l) as , where q 1 = k i − k k is the momentum transfer vector of the ERI.The 4h2p linear term calculation with the Madelung constant correction using a finite mesh K can be reformulated as using the change of variable k k → k i − q 1 and the periodicity of F kl (q 1 ).The error of this calculation compared to its TDL value, i.e., the first error term in Eq. (B9), can be written as Previously in [41], the quadrature error for F kl (q 1 ) with varying (k, l) is estimated as As to be demonstrated next, in the case of partially matched orbital indices (e.g., k = i, l = j), this error estimate turns out to be loose when k i = k a and can be further improved as where q = k a − k i + G 0 with G 0 ∈ L * chosen to make q ∈ Ω * .In the case of fully matched orbital indices (k = i, l = j), the Madelung constant correction is triggered in the ERI evaluation and can help remove the leading quadrature error when k i = k a as In the following discussion, we prove these two error estimates Eq. (B12) and Eq.(B13) at k i = k a separately.

Error estimate for linear terms with partially matched orbital indices
Consider a fixed set of (i, j, a, b) and k i , k j , k a ∈ K with k i = k a .Assume k = i and l = j.Our target is a sharper estimate of E Ω * F il (q 1 ), K q .Since K q and Ω * are both inversely symmetric over q 1 = 0, the quadrature error can be symmetrized as This symmetrized integrand can be further decomposed into two terms For the first term in Eq. (B14), we note that the ERI term can be detailed as The ERI nonsmoothness with q 1 ∈ Ω * comes from the first fraction term whose numerator is smooth and scales as O(|q 1 |) near q 1 = 0 using orbital orthogonality, It can be verified directly that H il eri (q 1 ) is singular at q 1 = 0 of order −1, and its symmetrization H il eri (q 1 ) + H il eri (−q 1 ) is singular at q 1 = 0 of order 0. Meanwhile, the amplitude H il amp (q 1 ) is smooth everywhere in Ω * except at q 1 = −q of order 0 according to Lemma 6.An error estimate lemma in [41] (restated as Lemma 24 in Appendix D) provides a quadrature error estimate for periodic functions in such a product form and can show that For the second term in Eq. (B14), direct application of Lemma 24 leads to error estimate of scale O(N ) is singular at q 1 = 0 of order −1.However, note that H il amp (q 1 ) is smooth at q 1 = 0 (recall that H il amp (q 1 ) is only singular at q 1 = −q = 0) and thus the subtraction H il amp (q 1 ) − H il amp (−q 1 ) scales as O(|q 1 |) near q 1 = 0. Multiplication by this extra O(|q 1 |) term improves the algebraic singularity of H il eri (q 1 ) at q 1 = 0 and the overall product can be shown to be singular at q 1 = 0 of order 0. Intuitively, this improved algebraic singularity at q 1 = 0 can lead to asymptotically smaller quadrature errors.To rigorously justify this statement, we generalize Lemma 24 to estimate the quadrature error for this special case.
x f 2 (0) = 0 for any derivative order |α| s.Assume γ > −d for f (x) to be integrable in V and γ + s + 1 0 so the leading algebraic singularity of f (x) is at x = 0. Consider an m d -sized uniform mesh X in V .Assume that X satisfies that z 1 , z 2 are either on the mesh or Θ(m −1 ) away from any mesh points, and m is sufficiently large that |z 1 − z 2 | = Ω(m −1 ).At x = z 1 and x = z 2 , f (x) is set to 0. The trapezoidal rule using X has quadrature error Remark 17.The two factors H d+1 V,z1 (f 1 ) and H d+1 V,z2 (f 2 ) characterize the algebraic singularities of the two functions, and their exact definition can be found in Appendix D. Proof of Lemma 16 is provided in Appendix D.2.
In order to utilize this result, we further decompose the second term in Eq. (B14) into H il eri (−q 1 ) H il amp (q 1 ) − H il amp (0) + H il eri (−q 1 ) H il amp (0) − H il amp (−q 1 ) .Applying Lemma 16 to both terms with γ = −1 and s = 0 gives where constant C is independent of i, j, a, b and k i , k j , k a ∈ K. Combining the above quadrature error estimates for the two terms in Eq. (B14), we prove a tighter error estimate at q = 0 shown in Eq. (B12) while the previous result in Eq. (B11) at q = 0 still holds.

Error estimate for linear terms with fully matched orbital indices
Let k = i, l = j and k i = k a and consider the corresponding calculation in the 4h2p linear term.The Madelung constant correction is triggered in the ERI evaluation at k k = k i or equivalently at q 1 = 0 and the calculation can be written as which uses the expansion of the Madelung constant with an arbitrary fixed parameter σ > 0 in Eq. (20).Note that the prefactor in the O(N −1 k ) remainder term above only depends on σ.The right hand side of the above reformation is equivalent to a singularity subtraction method that decomposes the original integrand into two terms, and then computes the numerical quadrature of the first term and the exact integral of the second term.The Madelung-corrected calculation thus has quadrature error only from the first term as Following a similar approach as in the partially matched case, the effective integrand after the Madelung constant correction above can be split into two terms as For the first term in Eq. (B15), the subtraction part can be detailed as where the numerator of the first fraction scales as O(|q 1 |) near q 1 = 0.This subtraction is thus singular at q 1 = 0 of order −1 and shares a similar form as the ERI with partially matched orbital indices.Using the earlier inverse-symmetry analysis for partially matched case, the quadarture error of the first term in Eq. (B15) when k i = k a can be estimated as For the second term in Eq. (B15), we exploit the inverse symmetry of K q and Ω * again and its quadrature error equals to that of its symmetrized version as This formula uses h σ (q 1 ) = h σ (−q 1 ) after symmetrization.Note that h σ (q 1 ) is singular at q 1 = 0 of order −2 while the term in the parenthesis scales as O(|q 1 | 2 ) near q 1 = 0 using the smoothness of H ij amp (q 1 ) at q 1 = 0. Therefore, the overall function above is singular at q 1 = 0 of order 0. To fit the integrand form in Lemma 16, we further decompose the symmetrized integrand above into Applying Lemma 16 to these two terms separately with γ = −2 and s = 1 gives Combining the above quadrature error estimates for the two terms in Eq. (B15), we prove a tighter error estimate at q = 0 shown in Eq. (B13) while the previous result in Eq. (B11) at q = 0 still holds, The same error estimate can be obtained for all the six linear terms that trigger the Madelung constant correction, and the remaining linear and quadratic terms contribute at most O(N −1 k ) quadrature error.Gathering the error estimates for all these terms together and plugging into Eq.(B9), we finish the proof.

B.3. Proof of Lemma 13: error in orbital energies
In the TDL, the orbital energy ε N k nk with any fixed n and k ∈ K converges to Comparing this exact orbital energy with its finite-size calculation in Eq. ( 2), the finite-size error without any corrections can be written as For the first quadrature error in Eq. (B16), i.e., the finite-size error in the direct term, the ERI with each i can be specified as Note that the momentum transfer vector of this ERI is always zero and the singular fraction term in this ERI is set to 0 by definition.As a result, the integrand is periodic and smooth with respect to k i ∈ Ω * .Therefore, the quadrature error of the direct term calculation thus decays super-algebraically according to Lemma 20 as For the second quadrature error in Eq. (B16), i.e., the finite-size error in the exchange term, the ERI with each i can be written as which is singular at q := k − k i = 0 of order −2 when n = i and 0 otherwise.An error estimate lemma in [41] (restated as Lemma 22 in Appendix D) gives a tight quadrature error estimate for such periodic functions with one point of algebraic singularity, and its application to the above integrand gives Combining the estimates of the two error terms in Eq. (B16), we obtain the overall finite-size error estimate for orbital energies without any corrections as From the above analysis, the dominant finite-size error in an occupied orbital energy lies in the calculation of the exchange term with i = n.In the Madelung-corrected orbital energy ε N k ,ξ nk , the correction is applied to this exchange term as (ignoring the prefactor −|Ω * | −1 ) Applying the change of variable k i → k − q and using the periodicity of the ERI with respect to k i , this corrected calculation can be reformulated as using the expansion of the Madelung constant with any fixed σ > 0 in Eq. ( 20).The quadrature error after the correction can be written as ).The effective integrand above after the correction can be detailed as The integrand singularity comes from the first fraction and is of order −1 since the ERI and h σ (q) share the same leading singular term.Similar to the analysis in the ERI contraction, we can then combine this singularity subtraction with the inverse symmetry of Ω * and K q to show that Combining this estimate with the above estimate for all the remaining direct and exchange terms in ε N k ,ξ nk , we obtain the finite-size error estimate for orbital energies with the Madelung constant correction as B.4. Proof of Lemma 14: error from previous iteration Fixing a set of entry index (i, j, a, b) and where the neglected terms are all the other linear and quadratic terms and the Madelung constant corrections to different terms are collected together at the end.In the subtraction A N k ,ξ (T ) − A N k ,ξ (S), the constant terms cancel each other.The subtraction between the two Madelung constant correction terms can be estimated directly as k ).The subtraction between the two 4h2p linear terms can be formulated and bounded as For the term with k k = k i (i.e., δ k k ,ki ), we have where the ERI definition at zero momentum transfer skips the singular fraction and is O(1), and according to the assumption on T − S it always holds that max ijab,ki ,kj,ka∈K For the term with k k = k a (i.e., δ k k ,ka ), its estimate is the same as the term above when k a = k i .When For the first summation term in Eq. (B17), we introduce the change of variable k k → k i − q 1 and write it as When q = 0, this term can be further bounded using Eq.(B18) by When q = 0, this term can be further bounded as where L * 0 denotes the set of 27 lattice vectors in L * around the origin.The first inequality is based on the lemma assumption on T − S that where G 0 ∈ L * 0 is the unique lattice vector that makes q 1 + q + G 0 ∈ Ω * .The third inequality can be obtained from the nonsmoothness characterization of function Using Lemma 11 in [41], f (z) is singular only at z = 0 of order −2 and its value at z ∈ Ω * \ {0} is bounded by C|z| −2 .
Based on the above estimates of the first term in Eq. (B17), we obtain the estimate of the error accumulation in the 4h2p linear term calculation as The same analysis can be applied to all the similar linear terms that contain ERIs with matched orbital indices.For other linear terms, the analysis can be done similarly and they all contribute at most O(N −1 k ) error to the overall subtraction.Taking the subtraction between the following 3h3p linear terms as an example, it can be formulated and bounded as For the subtraction between quadratic terms, we consider the 4h2p quadratic term as an example.The subtraction between the two 4h2p quadratic terms can be formulated and bounded as Similar analysis can be done to all the remaining quadratic terms, and they all contribute at most O(N −1 k ) error to the overall subtraction.Gathering all the estimates above together, we finish the proof.norm) instead, which is denoted by where T and t denote a generic amplitude computed using mesh K and in the TDL, respectively.Assume that the CCD amplitude equations using a sufficiently fine mesh K (with or without the Madelung constant correction) and in the TDL have unique solutions and denote the solutions as T N k * and t * respectively.Convergence of the CCD(n) amplitudes to the CCD amplitude is defined in the • 1 -norm sense as We impose a sufficient condition that guarantees the convergence of fixed point iterations by requiring the target mapping to be contractive in a domain that contains both the solution point and the initial guess.Assumption 18.For N k sufficiently large, the following statements hold: 1.The exact CCD(n) amplitude t n converges to the CCD amplitude t * point-wisely as n → ∞, i.e., and the initial guess 0, i.e., with a Lipschitz constant L < 1.
3. The exact CCD(n) amplitude t n and the domain B N k above satisfy that 4. For all the amplitudes {t n }, there exists a constant C such that Note that when consider the finite-size calculation with Madelung constant correction, the components ε N k or A N k in the second assumption need to be changed to ε N k ,ξ or A N k ,ξ accordingly.For the last assumption, we note that the two related error estimates in Lemma 11 have the prefactor C dependent on the amplitude t n .The assumption here is stronger in the sense that C needs to be independent of t n .k ).Proof.The finite-size error in the CCD energy calculation with or without the Madelung constant correction can be estimated as

Rigorous
where the last inequality uses the boundedness of |W ijab | in G N k and Lemma 10.To first estimate the above amplitude error when the Madelung constant correction is applied to both orbital energy and ERI contractions, we consider the error splitting Eq. (B2) for the amplitude calculation at the n-th fixed point iteration as , where the last estimate uses the assumption in Eq. (C4) for the first term and the assumptions of contraction maps in Eq. (C2) and M K t n−1 ∈ B N k in Eq. (C3) for the third term.Since the initial guesses in the finite and the TDL cases satisfy M K t 0 − T N k 0 1 = 0, we can recursively derive that and thus the first assumption Eq. (C1) gives Plugging this estimate into Eq.(C5) then finishes the proof for the scenario when the Madelung constant correction is applied to both orbital energies and ERI contractions.
For the two scenarios with partial Madelung constant correction, a similar analysis as above gives , where the dominant CN − 1 3 k error comes from uncorrected orbital energy or uncorrected ERI contraction term.Recursively, we can obtain the amplitude error estimate as which finishes the proof for the two scenarios with partial correction.
Appendix D: Quadrature error estimate for periodic functions with algebraic singularity This section presents a collection of lemmas that provide quadrature error estimates for trapezoidal rules over periodic functions with specific algebraic singularities, which are used in this paper.Most of the lemmas are proven in [41] and are restated here for completeness.In addition, we introduce and prove a new quadrature error estimate that is critical in describing the efficacy of the Madelung constant correction and inverse symmetry in reducing the quadrature error in orbital energies and ERI contractions.
The lemmas presented in this section provide the asymptotic scaling of the quadrature errors and also a quantitative relationship between the prefactors in the estimate and the algebraic singularities of the integrand.In addition to the singularity order, we further quantitatively characterize the algebraic singularity as follows.For a univariate function f (x) that is smooth everywhere in V except at x = x 0 with algebraic singularity of order γ, we define a constant .
For a multivariate function f (x, y) that is smooth everywhere in V X × V Y except at x = x 0 with algebraic singularity of order γ, we define a constant , where "•" in the subscript "(x 0 , •)" is a placeholder to indicate the smooth variable.Using these two quantities, we have following function estimates that will be extensively used in this section D.1.Existing results from Ref. [41] Lemma 20.Let f (x) be smooth and periodic in V = [− 1 2 , 1 2 ] d .The quadrature error of a trapezoidal rule using an m d -sized uniform mesh X in V decays super-algebraically as |E V (f, X )| C l m −l , ∀l > 0.
γ + s + 1 0, the three estimates together prove that the quadrature error of f (x) scales asymptotically as m −(d+γ+s+1) as m → ∞.In order to describe the extreme case where the two singular points are only O(m −1 ) away from each other for a given mesh X , i.e., δz = O(m −1 ), we provide a precise description of the three prefactors in the estimates above using δz.

Figure V. 1 .k 3 k
Figure V.1.Convergence of the CCD(1), CCD(2), CCD(3), and CCD correlation energies for a 3D periodic system of hydrogen dimers with increasing N k .The settings are distinguished by the presence or absence of the Madelung constant correction to the orbital energies ("ε N k ,ξ nk ") and the ERI contraction ("A N k ,ξ ").In (a)-(d), the dashed curves show the power-law fitting using C0 + C1N − 1 3

1 3k
example, if k i , k a ∈ K are adjacent to each other, |q| is of scale O(N − ) and the estimate in Lemma 11 suggests an error bound of scale O(N

Remark 19 .
The second assumption guarantees that {T N k n } lies in B N k and converges to T N k * .For the third assumption, Theorem 1 shows that with each fixed n > 0 the amplitudeT N k n converges to t n in the sense of lim N k →∞ M K t n − T N k n 1 = 0, suggesting that {T N k n } ⊂ B N k converges to t n with K → Ω * .Therefore B N k and M K t n should be related to each other, which leads to the third assumption.

1 3k
Statement of Theorem 2. Under Assumption 18, the finite-size error of the CCD correlation energy scales as O(N − ) in each of the following scenarios (1) the Madelung constant correction is only applied to the ERI contraction A N k (2) the Madelung constant correction is only applied to the orbital energy ε N k nk .When the Madelung constant correction is applied to both A N k and ε N k nk in the CCD calculation, the overall finite-size error scales as O(N −1 H l VX ×VY ,(x0,•) (f ) = min C : ∂ α x ∂ β y f (x, y) C|x − x 0 | γ−|α| , ∀|α|, |β| l, ∀x ∈ V X \ {x 0 }, y ∈ V Y = max |α| l ∂ α x ∂ β y f (x, y) /|x − x 0 | γ−|α| L ∞ (V ×V )
see Section IV.2.As a result, when applying the Madelung constant correction to both orbital energies and ERI contractions in the CCD(n) calculation, the overall finite-size error can be successfully reduced to O(N −1 k ).