Finite-size Effects in periodic EOM-CCSD for Ionization Energies and Electron Affinities: Convergence Rate and Extrapolation to the Thermodynamic Limit
Evgeny Moerman Alejandro Gallo Andreas Irmler Tobias Schäfer Felix Hummel Andreas Grüneis Matthias Scheffler


APS/123-QED




  
    Institute for Theoretical Physics, TU Wien,

    Wiedner Hauptstraße 8–10/136, 1040 Vienna, Austria 
  



  
    The NOMAD Laboratory at the FHI of the Max-Planck-Gesellschaft
  



moerman@fhi-berlin.mpg.de











§ ABSTRACT 
We investigate the convergence of quasi-particle energies
for periodic systems to the thermodynamic limit using increasingly 
large simulation cells corresponding to increasingly dense integration meshes in reciprocal space.
The quasi-particle energies are computed at the level of
equation-of-motion coupled-cluster theory for ionization
(IP-EOM-CC) and electron attachment processes (EA-EOM-CC).
By introducing an electronic correlation structure factor, the
expected asymptotic convergence rates for systems with different dimensionality
are formally derived.
We rigorously test these derivations through numerical simulations for 
trans-Polyacetylene using IP/EA-EOM-CCSD and the G_0W_0@HF approximation, 
which confirm the predicted convergence behavior.
Our findings provide a solid foundation for efficient schemes to correct 
finite-size errors in IP/EA-EOM-CCSD calculations.





  Finite-size Effects in periodic EOM-CCSD for Ionization Energies and Electron Affinities: 
  Convergence Rate and Extrapolation to the Thermodynamic Limit
    Matthias Scheffler 0000-0002-1280-9873
    September 9, 2024
=============================================================================================================================================================







§ INTRODUCTION


For most theoretical materials science studies, density functional
theory (DFT) is employed due to its favorable balance between
computational scaling and moderate accuracy. However, for many materials
properties the accurate inclusion of electronic exchange and correlation
effects is critical to achieve a qualitatively correct description for
scientifically and technologically important properties. In particular,
for the theoretical description
of electronic band gaps and band structures, most Kohn-Sham density functional
approximations (KS-DFAs) in use today are known to generally severely
underestimate
band gaps, oftentimes referred to as the band gap problem <cit.>.
Relaxing the condition of a multiplicative exchange-correlation potential of KS-DFT to allow
for a more flexible integral operator as it is the case in generalized KS-DFT resolves many of 
the issues associated with the band gap problem <cit.>, 
generally reducing the band gap error and rectifying the missing derivative discontinuity.



Apart from hybrid functionals, which are the most widely used representatives 
of generalized KS-DFAs, the electronic structure method of choice for the
calculation of quasi-particle energies has become the
GW-approximation <cit.>, which takes the DFA
one-electron wave functions as a starting point but explicitly
accounts for electronic exchange and
correlation effects using a perturbation theory approach.
The GW-approximation yields significantly improved band gaps
compared to the most widely used approximate density functionals.
However, the most commonly used method based on the GW approximation, the
G_0W_0 method,
is known to have its limitations as well, the most glaring one being
the dependency of the band gap result on the underlying DFA, the so-called
starting point dependence <cit.>.
Other higher-order corrections to the GW approximation require
the inclusion of so-called vertex corrections in the self-energy and
the screened interaction W <cit.>.
Although certain improvements can be achieved
using vertex corrections, it remains challenging to systematically improve
upon the GW approximation, which already achieves an excellent trade-off
between computational cost and accuracy <cit.>.

A systematically improvable method is the equation-of-motion
coupled-cluster <cit.>(EOM-CC) framework.
EOM-CC theory, being the extension of ground-state CC theory to excited states, allows to
theoretically describe systems upon removal (IP-EOM-CC), addition
(EA-EOM-CC) or vertical excitation (EE-EOM-CC) of an electron. As the
electronic band gap is defined as the difference between the electron
affinity/attachment (EA) and the ionization potential (IP), it is possible
to obtain band gaps and entire band structures in the EOM-CC
framework <cit.>.
The relation between EOM-CC and GW band gaps was also investigated,
showing the differences and similiarites between the diagramatic contributions 
and the results of both approaches <cit.>.
However, the major obstacle of ab initio calculations employing EOM-CC theory,
and CC theory in general, is the high computational cost and excessive
memory requirements. In addition to that, CC theory
explicitly incorporates long-range electronic exchange and correlation contributions,
making convergence to the bulk-limit significantly slower than it is the case for DFAs.

Due to the
high computational cost of CC methods, it is even more challenging
than for DFA and GW calculations to converge to the tdl, 
which is approached by increasing the number of particles N_part in the 
simulation cell, N_part→∞,  while keeping the particle density constant. 
The GW approximation partly ameliorates this problem by adding corrections
for the long-range behavior of the dielectric function using
k· p-perturbation theory, which are often referred to as
head- and wing-corrections <cit.>. 
In recent years,
studies of electronic band gaps via Quantum Monte Carlo (QMC)
methods have been published as well <cit.>,
where finite-size effects were also discussed as one of the major sources of
error. For these QMC band gaps, the N-electron ground-state
and the electronic state with one electron more (N+1) or less (N-1) were
determined separately and the energy difference of these states was computed
to obtain the band gap value. As the leading-order contribution to the
finite-size error, the interaction of the added particle (N+1) or
hole (N-1) with its periodic images was identified, which was corrected
by subtracting the screened Madelung term from the quasi-particle band gap.
Higher-order finite-size errors resulting from multipole moments of the
charged states were corrected by means of system size extrapolation.
Unfortunately, these approaches are not straight-forwardly applicable to
EOM-CC methods: The dielectric function, and therewith the head- and
wing-correction to it, is not directly accessible in the canonical
EOM-CCSD formulation. For that, linear response CC theory would be
necessary <cit.>. In contrast to the QMC ansatz,
CC methods do not use trial wavefunctions but generally rely on the
hf wave function as a starting point, which already
incorporates the Madelung term. Instead, one needs to resort to
extrapolation techniques
or more sophisticated
finite-size error estimations using the transition structure
factor <cit.>, which,
however, has been only
formulated for ground-state CC theory so far and is only a viable
technique if the correlation structure factor is represented in a 
plane wave (PW) basis.
Using a very different approach for the simulation of crystalline systems
in the CC framework, it has been
demonstrated that via  a cluster embedding approach accurate band gap
predictions can be achieved <cit.>. It must be stressed,
however, that the results discussed in the present work
assume periodic boundary conditions. A second point of departure is that
the work on band gaps from cluster embedding techniques utilized similarity transformed
EOM (STEOM) theory, while we in this work explore the EOM-CC method.

In this work, all CC and EOM-CC calculations were performed using a super cell approach.
Even though
it is entirely sufficient and – due to the exploitation of the translation symmetry –
computationally significantly more efficient to calculate the band gap of a perfect crystal
on a regular k-grid in reciprocal space using the primitive unit cell, a k-point aware,
block-sparse treatment of the CC and EOM-CC equations is not yet available in CC4S.
Even for DFAs, it is well known that the super cell size convergence is an impractical 
approach for solving a proper k-point summation.
The IP- and EA-EOM-CCSD method exhibit a computational scaling of
N_o^3N_v^2N_k^3 and N_oN_v^4N_k^3 <cit.>, respectively, where N_o, N_v and N_k denote
the number of occupied and unoccupied orbitals per unit cell and the number of k-points.
The memory scaling is dominated by the T_2-amplitudes and is proportional to
N_o^2N_v^2N_k^3. If instead of N_k k-points, a super cell approach with
N_u=N_k unit cells is employed, a (N_uN_o)^3(N_uN_v)^2 = N_o^3N_v^2N_u^5
and a (N_uN_o)(N_uN_v)^4 = N_oN_v^4N_u^5 computational scaling for the
IP- and EA-EOM-CCSD method, respectively, is the consequence. Analogously,
the memory scaling becomes proportional to N_o^2N_v^2N_u^4.
Thus, a k-point based treatment of the EOM-CC equations
would result in a reduced computational scaling of a factor
N_k^2 and a reduction in memory scaling by a factor
of N_k. 
It must, however, be stressed that in accordance with
Bloch's theorem, the bvk cell of a M× K× L
super cell evaluated at a single k-point k_off
is identical to the bvk cell resulting from a
primitive unit cell being evaluated on a regular M× K× L
k-grid shifted by k_off. Hence, even though the super cell
approach is notably more computationally expensive, the numerical EOM-CCSD
results presented here are not affected by this choice. Still, as is well-known
for standard electronic-structure theory, only the k-summation approach
is practically feasible in order to achieve covergence. 

For the EOM-CC methods, one is currently forced to perform
calculations of increasing system size and perform an extrapolation
to the tdl <cit.>, which requires knowledge
about the convergence rate of the correlation energy with respect to
the system size.  While this convergence rate has been studied for the ground-state of the
3-dimensional case of a bulk solid <cit.>,
the formal convergence behavior for electronically excited states (charged or neutral) for any dimension is unknown.


In this work, we formally derive the analytical expression governing the convergence
rate of the band gap on the IP/EA-EOM-CCSD level of
theory. Subsequently, we verify the correctness of the derived
expression by applying it to the band gap of a single chain of
trans-Polyacetylene, demonstrating an efficient 
extrapolation approach to the tdl. 
Furthermore, by repeating the
calculations using
the G_0W_0 method with a hf starting point (G_0W_0@HF), we show that
the derived convergence rate is not specific to periodic EOM-CC theory
but can be used for other correlated methods as well.






§ THEORY



 §.§ The EOM-CC theory

The EOM-CC framework is an extension to ground-state CC
theory, to compute properties of excited states.  Depending on the
nature of the excitations, different EOM-CC methods are available: The
most prominent ones are EE-EOM-CC, for neutral electronic excitations,
IP-EOM-CC for ionization processes and EA-EOM-CC for electron
attachment processes. For the present work only the latter two methods
are of importance, as the fundamental band gap of a material is given
by the difference of its lowest ionization potential and electron affinity.

The starting point of the EOM-CC method is the ground-state CC
many-electron wave function |Ψ_0⟩, which is defined by an
exponential ansatz


    |Ψ_0⟩ = e^T̂|Φ_0⟩,


with the Slater determinant |Φ_0⟩, which is the ground-state wave
function of a preceding mean-field calculation, usually hf.
T̂ = ∑_n^MT̂_n is the so-called cluster
operator, which can excite up to M electrons:


    T̂      =T̂_1 + T̂_2 + ⋯+T̂_M
          =∑_i,at^a_iâ^†_aâ_i + ∑_i,j,a,b1/4t^ab_ijâ^†_aâ^†_bâ_jâ_i + ⋯
          +(1/M!)^2
        ∑_i,j⋯
     a,b⋯
        t^ab⋯_ij⋯â^†_aâ^†_b⋯â_jâ_i


with the coefficients t^a_i,
t^ab_ij,⋯ being again the
cluster amplitudes and â^†_p and â_q
the creation/annihilation operators
in second quantization, creating/annihilating an electron in orbital p/q.
The notation in Equation <ref> is
such that i,j,k
denote occupied  and
a,b,c unoccupied spin orbitals. 
If M is equal to the number of electrons N of the system, the
ansatz in Equation <ref> is exact. For reasons
of otherwise impractical computational scaling, T̂ is truncated
in practice.
For example, M=2 corresponds to CC theory
with single- and double excitations (CCSD),
which is also used in this work.

To compute the wave function of the n-th
excited state, the EOM-CC framework
starts from a linear ansatz


    |Ψ_n⟩ = R̂_n|Ψ_0⟩


where R̂_n is an excitation operator
similar to the cluster operator T̂.
For IP-EOM-CC and EA-EOM-CC, R̂_n
assumes the form


    
    R̂^IP_n =
            ∑_i r_i,nâ_i +
            ∑_ija r^a_ij,nâ^†_aâ_jâ_i + ⋯

    
    R̂^EA_n =
            ∑_a r^a_nâ^†_a +
            ∑_iab r^ab_i,nâ^†_aâ^†_bâ_i + ⋯.



The r-coefficients contain the description of the wave function of
the n-th excited state and will be summarized under the excitation vectors
|R_n^IP⟩ or |R_n^EA⟩ .
The objective of the EOM-CC methodology
is to determine these coefficients.  In analogy
to the T̂-operator, the operators defined in Equations <ref>
and <ref> contain all possible
processes to excite an N-electron system into an (N-1)-state and
an (N+1)-state, respectively. However, for the same reasons stated
before for ground-state CC, R̂^IP_n and
R̂^EA_n are truncated in practice. The most common
approximation for EOM-CC involves only the first two terms of Equations
<ref> and <ref> and is termed
IP-EOM-CCSD and EA-EOM-CCSD, respectively.  Consequently, IP-EOM-CCSD
accounts only for 1-hole- and 2-hole-1-particle-excitation processes,
while EA-EOM-CCSD is restricted to 1-particle- and 2-particle-1-hole
processes. We stress that, here, the term excitation process is not restricted
to charge neutral processes but includes the removal and addition of
electrons.

In order to determine the excited states |R_n^IP/EA⟩ and
the related IP- or EA-energy, the eigenproblem


   
    H̅|R_n^IP⟩ =
            IP_n|R_n^IP⟩

    
    H̅|R_n^EA⟩ =
            EA_n|R_n^EA⟩



needs to be solved, where H̅ = e^-T̂Ĥe^T̂
is the similarity-transformed Hamiltonian, with 
T̂ being the ground-state CC cluster operator.
Due to the
generally intractable size of H̅ in the chosen representations,
the eigenvalues and -vectors cannot be determined directly but need to be computed via an
indirect approach like Davidson's method <cit.>.

Once the excited states |R_n^IP/EA⟩ are obtained, the
corresponding IPs or EAs can be computed via


    IP_n/EA_n =
        ⟨ R_n^IP/EA|H̅|R_n^IP/EA⟩/⟨ R_n^IP/EA|R_n^IP/EA⟩.


Note here, that even though
H̅ is for the given T̂ non-symmetric and therefore every eigenvalue
is associated with both a left and right eigenvector,
the calculation of the eigenvalue in
Eq. (<ref>) only requires the knowledge of one of the eigenvectors.
The working equations for H̅|R_n^IP/EA⟩ can be
found in, e.g. <cit.>.




Finally, note that while the theoretical framework of CC theory for the ground
and excited states was elucidated using spin orbitals, henceforth all quantities
will be expressed in terms of (spin-independent) spatial orbitals, as all
the results shown in this work have been obtained without the consideration of
spin degrees of freedom.



 §.§ The EOM-CC structure factor

Let us now introduce an expression for the IPs and EAs
that makes it possible to analyze their dependence on the interelectronic
distance.
A closer look of the working
equations show
that all contributions to the expectation value in
Eq. (<ref>) constitute contractions of cluster
amplitudes (see Eq. (<ref>)), r-coefficients
(Eqs. (<ref>)/(<ref>)) and Coulomb
integrals


    V^pq_rs = ∬d𝐫d𝐫'
        ϕ_p^*(𝐫)ϕ^*_q(𝐫')
        ϕ_r(𝐫)ϕ_s(𝐫')/|𝐫-𝐫'|,


where ϕ_p(𝐫) denotes a single-particle
state of the underlying mean-field
theory (hf in this work).

By assuming bvk boundary conditions we can introduce the co-densities in reciprocal
space

    C^p_r(𝐪) =
        ∫d𝐫ϕ^*_p(𝐫)ϕ_r(𝐫)e^-i𝐪𝐫.


One can rewrite Eq. (<ref>) as


    V^pq_rs =∑_𝐪 w_𝐪C^r_p^*(𝐪)v(𝐪)C^q_s(𝐪),


The discrete 𝐪-vectors in Equation <ref> and <ref>
lie on a grid in reciprocal space and
v(𝐪) = 4π/|𝐪|^2 is
the Coulomb potential in reciprocal space for the three-dimensional case.
We stress that the 𝐪-mesh is used to represent the co-densities in Fourier space.
If the single-particle states are expressed using Bloch's theorem such that ϕ_s(𝐫)=e^i𝐤_s𝐫u_s(𝐫),
where u_s(𝐫) is a cell periodic function and 𝐤_s is a wave vector in the first Brillouin zone,
the 𝐪-vectors correspond to the difference between the corresponding
wave vectors and a reciprocal lattice vector of the periodic unit cell.
If the single-particle states are represented using a supercell approach, the 𝐪-vectors correspond to
reciprocal lattice vectors of the supercell. Both approaches are formally equivalent, although
Bloch's theorem enables a computationally more efficient implementation.
w_𝐪 is a weighting factor that depends on the employed integration grid and method.
It follows from Eq. (<ref>)
that C^q_s(0)  is equivalent to the
overlap between the two involved single particle states

    C^p_r(𝐪=0) = δ_p,r.



The terms contributing to the expectation value of
IP- and EA-EOM-CCSD of Equation <ref>, that is the IP and EA, can be broadly
separated into two types, single-body and many-body
contributions.

We define single-body mean-field contributions to explicitly depend on the
Fock matrix elements f^p_q and only contain Coulomb integrals
implicitly by virtue of the Hartree and exchange contribution of the Fock matrix elements. A
representative single-body mean-field contribution to the
IP and EA is given by



    
    IP_n = ⋯ -r^ij*_b,n f^l_jr^b_il,n+⋯

    and
    
    EA_n = ⋯ +r^j*_ab,n f^a_c r^cb_j,n+⋯.



We define IP^(1)_n/EA^(1)_n as the sum of all single-body mean-field terms
in the expression for IP_n/EA_n.
It should be noted that the commutator expansion of the similarity transformed Hamiltonian
also gives rise to effective single-body contributions originating from the contraction
of the two-body Coulomb operator with different orders of T̂_1 and T̂_2.
However, in this work we choose
to include only terms from the underlying mean-field Hamiltonian in our definition
of single-body terms.

It follows from our definition of IP^(1)_n/EA^(1)_n that all remaining
contributions to IP_n/EA_n are referred to as many-body terms and
depend explicitly on Coulomb integrals given by, for example,


    
    IP_n = ⋯
            -2r^i*_n V^kl_cd t^cd_il r_k,n
            +⋯

    and
    
    EA_n = ⋯
            -2r^*_a,n V^kl_cd t^ad_kl r^c_n
            +⋯,



where r^i_n,r_a,n and r^a_ij,n, r^ab_i,n
denote the single- and double excitation coefficients
of the n-th IP- and EA-EOM-CCSD excitation operator |R_n^IP/EA⟩, respectively.
A diagrammatic representation of two exemplary many-body EA-EOM-CCSD contributions is shown in Figure <ref>.

By replacing all explicit occurrences of the Coulomb integrals V^pq_rs in the many-body contributions
 by the decomposition in Equation (<ref>)
and by contracting over all particle- and hole-indices
involved in the evaluation of Eq. (<ref>),
one arrives at an expression for
the IP_n and EA_n expectation value as a
sum of IP^(1)_n/EA^(1)_n and a product
of the EOM-CC structure factor S_n^IP/EA(𝐪) with the Coulomb potential v(𝐪) as
shown in Eqs. (<ref>) and (<ref>)


  
    IP_n = ∑_𝐪w_𝐪 S^IP_n(𝐪)v(𝐪)+IP^(1)_n

  
    EA_n = ∑_𝐪w_𝐪 S^EA_n(𝐪)v(𝐪)+EA^(1)_n,



where the first term of Equation <ref> (Equation <ref>) contains all the many-body,
or “correlation” contributions to the IP_n (EA_n) as introduced in Equation <ref>,
while the second term IP^(1)_n (EA^(1)_n) denotes the sum of all single-body
mean-field contributions defined in Equation <ref>.

The above expression gives access to the dependence of the “correlation” energy
contribution to IP_n/EA_n on the momentum transfer vector 𝐪,
making it possible to perform a Fourier transform of S^IP/EA_n(𝐪)
and attribute contributions to IP/EA_n
on the inter-electronic distance.
In this work, however, we restrict ourselves to studying the dependency on
the momentum transfer vector 𝐪.







 §.§ The EOM-CC structure factor at 𝐪=0







Since C^q_s(0)=δ_q,s, the contributions to
S^IP/EA_n(𝐪=0) can be further simplified.
In particular, one discovers that
S^IP/EA_n(𝐪=0) corresponds to the many-body character of the
IP or EA state. We define
the single-body character p_1 of an IP or EA state as


    
    p_1,n = ∑_i |r_i,n|^2

    and
    
    p_1,n = ∑_a |r^a_n|^2,



respectively, where the r-coefficients
from Equation <ref> and <ref> are
used. In this definition, the single-body character p_1,n of the
n-th IP or EA state is given by the contribution of single-hole or
single-particle processes to the overall description of the (N+1)/(N-1) wave
function. For a normalized excited state |R_n^IP/EA⟩,
p_1,n lies between 0 and 1.
As we will show, for both IP- and EA-EOM-CCSD,

    S^IP/EA_n(𝐪=0) =
        -(1-p_1,n)

holds. Thus, the value of the
EOM-CC structure factor at 𝐪=0 is equal to the negative
many-body character of the excitation.

This results from Equation <ref>, so that the only terms of the EOM-CCSD
equations, that contribute to the value
of the EOM-CC structure factor at 𝐪=0 are those
which feature Coulomb integrals of type V^ab_cd,
V^ij_kl, V^ai_bj or V^ia_jb. That is,
the integrals must be representable as products of hole-hole or particle-particle co-densities following Equation
<ref>. In the case of IP-EOM-CCSD,
there are 5 terms which contribute to S(𝐪=0), which are


    S^IP_n(𝐪=0)=
         [  .-   r^ij*_b C^l*_j(𝐪) C^b_d(𝐪)r^d_il
    
        -   r^ij*_b C^k*_i(𝐪) C^b_d(𝐪)r^d_kj
       r^ij*_b C^k*_i(𝐪) C^l_j(𝐪)r^b_kl
    
        -   r^ij*_b C^k*_i(𝐪)C^b_d(𝐪) t^d_j r_k
       r^ij*_b C^k*_i(𝐪) C^l_j(𝐪) t^b_l r_k. ]_𝐪=0.


By making use of the fact that at 𝐪=0, the co-densities reduce to overlap integrals between the single-particle states (see Equation <ref>),
Equation <ref> reduces to


    S^IP_n(𝐪=0)=-   r^ij*_br^b_ij
        -r^ij*_b r^b_ij
        +r^ij*_b r^b_ij
    
        -   r^ij*_bt^b_j r_i
        +r^ij*_b t^b_j r_i
    
        =-   r^ij*_br^b_ij,


which – for a normalized EOM-CCSD excitation vector –
is equivalent to Equation <ref>. The derivation for the EA-EOM-CCSD structure factor is analogous.

We emphasize that the value for S^IP/EA_n(𝐪=0) is a direct
consequence of the chosen definition of IP^(1)_n/EA^(1)_n
and the eigenvalue in Eq. <ref>. Although this choice is
helpful for the analysis of finite-size errors, there also exist
alternative approaches, for example, using left eigenvectors as bra-states in
Eq. <ref>, resulting in a different value for p_1,n
and S^IP/EA_n(𝐪=0).
We expect that for the systems studied in this work, which exhibit a very small
many-body character, both approaches discussed above would yield very similar results.



 §.§ Long-range behavior of the EOM-CC structure factor

To determine the asymptotic behavior of the
EOM-CC structure factor in the long-wavelength limit, that is for |𝐪|→ 0,
one can perform a Taylor series expansion around 𝐪=0 by computing the
derivatives
of the EOM-CC structure factor with respect to 𝐪.
The only explicit 𝐪-dependence of that quantity comes from
the co-densities C^p_r(𝐪).
As laid out in Section
<ref>, v(𝐪) is not included in the expression for
S^IP/EA(𝐪), so that the products of co-densities
C^p*_r(𝐪)C^q_s(𝐪) are the only
𝐪-dependent quantities. We note that
we neglect the implicit 𝐪 dependence possibly introduced by the EOM-CC amplitudes
(r_i, r^a, r_ij^a, r_i^ab). This assumed 𝐪 independence of the 
r-amplitudes for 𝐪→ 0 is justified by the fact that for systems with a band gap none of the quantities 
that appear in the EOM-CC equations, which are the Fock matrix elements, the Coulomb integrals 
and the ground-state T-amplitudes, depend on 𝐪 (assuming the crystal momentum is conserved).
While this is evident for the Fock matrix and the Coulomb integrals, the 𝐪 independence
of the T_2-amplitudes is necessary for the ground-state transition 
structure factor to yield the correct 1/N_k convergence behavior
in the thermodynamic limit (𝐪→ 0) <cit.>. For the T_1 amplitudes,
which do not explicitely appear in the expression of the ground-state CC correlation energy or the
transition structure factor, a similiar argument as for the EOM-CC equations applies: Since the
T_1-amplitudes are determined iteratively from 𝐪 independent quantities, it is 
reasonable to assume that they themselves do not exhibit any dependence on 𝐪 either. 

By making use of the product rule and by realizing
that a co-density at 𝐪=0 is equivalent to the
overlap between the two involved single particle states as defined by Eq. (<ref>),
one finds that


    ∂/∂𝐪 C^p*_r(𝐪)C^q_s(𝐪)|_𝐪=0 =
        δ_p,r∂/∂𝐪 C^q_s(𝐪)|_𝐪=0 +
    δ_q,s∂/∂𝐪 C^p*_r(𝐪)|_𝐪=0


Eq. (<ref>) reveals, that
if the co-densities and by extension the Coulomb integrals in the
EOM-CCSD equations would only couple holes with particles, the first
derivative of the EOM-CC structure factor would vanish.
In passing we note that this is the situation
for the ground-state CC transition structure factor, which is why the
lowest 𝐪-order in the long-wavelength limit for insulators is
quadratic <cit.>.  The EOM-CC working equations,
however, do also feature contractions with Coulomb integrals, which
couple holes with holes and particles with particles, so that formally
a linear contribution in 𝐪 to the EOM-CC structure factor must be considered,
because in general
∂/∂𝐪 C^q_s(𝐪)|_𝐪=0≠ 0.

By repeating this procedure for the second derivative,
one finds that the quadratic contribution
does in general not vanish either, so that we can approximate
the asymptotic behavior of the EOM-CC
structure factor up to second order


    lim_|𝐪|→ 0  S^IP/EA(|𝐪|) ≈
        S^IP/EA_n(𝐪=0) + α·|𝐪| + β·|𝐪|^2,


where α and β are constants and S^IP/EA_n(𝐪=0) is the value
of the EOM-CC structure factor at 𝐪=0 as discussed in Section <ref>.
We note that in Equation <ref>,
we assume S^IP/EA_n(𝐪) to be spherically symmetric, which – while not generally the case – simplifies the 
following derivation and analysis.

Note, that the above derivation – while based on some assumptions– did not 
depend in any way on the dimensionality of the system.




 §.§ Convergence rate of IPs, EAs and band gaps
to the TDL


Based on the discussion from the previous sections, we now want to make a statement about the asymptotic
convergence behavior of computed IPs, EAs and band gaps to the tdl.
The tdl is approached for supercells with increasing size or increasingly dense k-mesh used to sample
the first Brillouin zone.
This corresponds to a continuous  𝐪-representation
of the EOM-CC structure factor such that Eqs. (<ref>) and
(<ref>) are given by


    IP_n/EA_n = ∫d𝐪  S^IP/EA_n(𝐪)v(𝐪) + IP_n^(1)/EA_n^(1),


Finite size errors of, for example, IP/EA energies are defined by the difference between calculations
using finite system sizes and the tdl.
For the above equation this corresponds to the difference between a continuous integration and a
discrete sampling. Similar to the case of ground state CC calculations, we assume that the
values of the EOM-CC structure factor at the sampled 𝐪 points converge rapidly with 
respect to the employed system size, which will be later verified numerically.
In other words, S^IP/EA_n(𝐪)=S^IP/EA-(TDL)_n(𝐪),
where S^IP/EA_n(𝐪) and S^IP/EA-(TDL)_n(𝐪) refer to the structure
factor from the finite system and tdl, respectively. Note that in this approach 𝐪 is restricted to a
discrete subset depending on  the system size, whereas its continuous in the tdl.
In this case, the largest contribution to the finite size error originates from the
employed finite simulation cell size and the neglect of long-range interactions in real space
corresponding to short 𝐪-vectors.


Under these assumptions it is reasonable to define the following estimate for the finite size
error of the correlation contribution Δ_FS^IP/EA


    Δ_FS^IP/EA =
        ∫_Ω_q_mind𝐪 S^IP/EA_n(𝐪)v(𝐪).


The integral in Equation <ref> is carried out over
a sphere Ω_q_min centered at the Γ point with radius q_min. The radius
should be understood as a measure for the shortest reciprocal lattice vector of
the considered simulation cell.
For the 3-dimensional case, we evaluate the integral
in Equation <ref> using the Fourier
transform of the Coulomb potential in three dimensions which is


    v(𝐪) = 4π/𝐪^2,


so that – using q=|𝐪| – a measure of the finite-size error is given by


    Δ_FS^IP/EA   ∝∫_Ω_q_mind𝐪 S^IP/EA_n(𝐪)v(𝐪)
       ∝∫_0^q_mind q  4π q^2
        S^IP/EA_n(𝐪)v(𝐪) 
       ∝∫_0^q_mindq  S^IP/EA_n(𝐪=0) + α· |𝐪|+ β· |𝐪|^2 
       ∝ 
        A· q_min +
        B· q^2_min +
        C· q^3_min,


where all constants resulting from the integration
are collected in the parameters A, B and C. By assuming spherical symmetry of the
EOM-CC structure factor, the derivation in Equation <ref>
is simplified substantially. 















More practical is Equation <ref> if expressed as a function
of the total number of 𝐤-points N_k of a uniform k-mesh.
Noting, that


    q_min∝ N_k^-1/d,


where d is the dimension of the system. By applying this
equality for the three-dimensional case, Equation <ref> becomes


    Δ_FS^IP/EA   ∝
        A· N_k^-1/3 +
        B· N_k^-2/3 +
        C· N_k^-1.




  §.§.§ Convergence rates for low dimensional systems


Since we assume that the asymptotic behavior of the EOM-CC structure factor is
independent of the dimension of the system, we can now make use of
Equation <ref> to derive the leading-order
rate of convergence to the tdl for one- and two-dimensional systems.
To this end we make use of the
Fourier transforms of the Coulomb potential in the two- and one-dimensional case <cit.>:


    v^2D = 2π/|𝐪|


and


    v^1D = -2γ_E + ln(1/𝐪^2),


respectively, where γ_E denotes the Euler constant. 
Note that we employ here the lower-dimensional Fourier transforms
of the Coulomb potential even though realistic one-dimensional systems 
generally have a finite extent in the directions perpendicular to
the a one- or two-dimensional material.
In the long-wavelength limit, however, which is the focus of the present work,
these contributions become negligible in comparison to those of the
sheet or chain direction(s). This justifies the representation
of the Coulomb potential by its lower-dimensional Fourier transforms.
By repeating the same steps as for 
the derivation of the 3-dimensional convergence rate, 
we find that the finite-size error in a two-dimensional
system converges like


    Δ_FS^IP/EA   ∝
        A· N_k^-1/2 +
        B· N_k^-1 +
        C· N_k^-3/2.


and for a one-dimensional system the convergence rate is
given by


    Δ_FS^IP/EA(N_k)    ∝    A· N_k^-1ln(N_k^-1) + B· N_k^-1
        +     C· N_k^-2ln(N_k^-1) + D· N_k^-2
        +     E· N_k^3ln(N_k^-1) +
        F· N_k^-3.


A summary of the derived convergence rates for
1, 2 and 3 dimensions is given in Table <ref>.
We reiterate that these convergence rates are estimated from the contributions
to IP_n/EA_n around 𝐪=0 in a sphere
with a radius decreasing as N_k increases.
It should also be noted that the actual convergence rate of numerically
computed IP_n/EA_n's depends on the chosen
treatment of the Coulomb singularity in the respective computer
implementation, which will be discussed in the following section.























 §.§ Convergence rates and treatment of Coulomb singularity


The convergence rates derived in the previous section
assume a spherically truncated integration around 𝐪=0 to estimate
the finite size errors Δ_FS^IP/EA.
In practical ab initio calculations, however, there exist a variety
of treatments to approximate the integral around the Coulomb singularity at 𝐪=0,
which strongly influence the convergence rates as can already be observed for the exchange energy 
contribution <cit.>.
We now briefly discuss the significance of  the treatment of the Coulomb singularity for
the convergence rates derived in the previous section.

For the present study, we employ a Coulomb singularity treatment
that captures the contribution of the
S^IP/EA_n(𝐪=0) term to the
integral in Eq. <ref> exactly if
the value of S^IP/EA_n(𝐪=0) is converged with respect to system size.
As a consequence, the contributions to the finite size errors proportionate to A given
in Table <ref> will already be accounted for and the expected
next-leading order contribution to the  finite size error will be proportionate to B.
In particular, the plane wave basis set calculations of the present work compute
the average Coulomb kernel for the volume element at 𝐪=0 to estimate
its contribution to the integral <cit.>.

We note that there also exist other approaches in the literature that disregard the
Coulomb singularity contribution



and obtain EOM-CC band gaps by extrapolation.
Refs. <cit.>
disregard the 𝐪=0 contribution already in the underlying hf calculation.
As a consequence the hf band gap is underestimated
and converges only as N_k^-1/3 in three-dimensional systems.
The same applies to the finite-size error in EOM-CC theory as shown in the previous section.
It should be noted that the finite size errors from hf and EOM-CC partly cancel each other.
However, the extrapolation procedure still requires a careful checking and can lead to errors
that are difficult to control due to next-leading order contributions to the finite size errors.














The last piece of information required to
compute this integral is the
Fourier transform of the Coulomb potential 1/𝐫 in 1 dimension, which is


    ∫ d𝐫 e^i𝐪𝐫/𝐫 = -2γ_E + ln(1/𝐪^2),


γ_E being the Euler constant.

Inserting Equation <ref> and
Equation <ref> into Equation
<ref>, yields


    Δ_FS^IP/EA   =
        ∫_0^q_min
        d𝐪 S^IP/EA_n(𝐪)v(𝐪)
       ∝∫_0^𝐪_min
        dq [S(𝐪=0) + a· |𝐪|+b· |𝐪|^2][-2γ_E + ln(1/𝐪^2)]
       ∝ A· q_minln(q_min) +
        B· q_min + 
       C· q_min^2ln(q_min) +
        D· q_min^2 + 
       E· q_min^3ln(q_min) +
        F· q_min^3,


where all constants resulting from the integration
are collected in the parameters A, B, C, D, E and F.

More practical is Equation <ref> if expressed as a function
of the number of 𝐤-points N_k along
the periodic dimension. Noting, that


    q_min∝1/N_k


Equation <ref> becomes


    Δ_FS^IP/EA(N_k) ∝   A·1/N_kln(1/N_k) +
        B·1/N_k + 
       C·1/N_k^2ln(1/N_k) +
        D·1/N_k^2 + 
       E·1/N_k^3ln(1/N_k) +
        F·1/N_k^3.







 §.§ Finite-size convergence for 1, 2 and 3 dimensions

Since the asymptotic behavior of the EOM-CC structure factor is
independent of the dimension of the system, we can make use of
Equation <ref> to derive the leading-order
rate of convergence to the tdl for 1, 2 and 3 dimensional systems.
By following the same steps as in Equation Section <ref>,
we obtain these convergence rates, which are summarized in Table <ref>.







§ COMPUTATIONAL DETAILS

As the practical representation of the EOM-CC structure factor requires the utilization
of a plane wave basis, we employed the
Vienna Ab Initio Simulation Package (VASP) <cit.>
in combination with the CC4S software package, a periodic CC
code <cit.>, where the working equations of the IP- and
EA-EOM-CCSD methods and their respective structure factors were implemented.
In the case of the LiH EOM-CC and ground-state CC structure factors, which are discussed in 
Section <ref>, a relatively small basis set with 3 virtual orbitals per occupied 
orbital was employed in order to compute the huge super cell sizes necessary to resolve 
the structure factor for very small momentum vectors.
However, we are concerned with the qualitative long-range description of the
EOM-CCSD structure factor, for which such a small plane-wave basis set is sufficient.
Similarly, the EOM-CC structure factors
for trans-Polyacetylene (tPA) were computed using a reduced basis set of 4 virtual orbitals
per occupied orbital. Due to substantial amount of vacuum in the simulation cell of the tPA chain,
it is not practical to converge the results with respect to the plane-wave energy cut-off, which is why
VASP was only used to obtain the EOM-CC structure factor but not the electronic band gaps themselves.
Specifically for the LiH calculations, the  and  POTCARs were employed, using
a plane-wave energy cut-off of 300 eV. The tPA EOM-CCSD structure factor was computed with the  and  POTCARs and an energy cut-off of 300 eV.

To alleviate convergence
problems of the hf and the CC calculations associated with
the discretization of the Brillouin zone (BZ) for the highly anisotropic simulation cells of tPA, we employ the recently
developed improved sampling method of the Coulomb potential
in VASP <cit.>.


The numerical convergence tests for the band gaps of tPA
were performed using FHI-aims <cit.>,
which employs numeric atomic orbitals and which has been interfaced to CC4S as well <cit.>.
To determine the convergence to the tdl, the unit cell of a single tPA chain with 2 carbon atoms and 2 hydrogen atoms
with at least 80 A of vacuum in each direction perpendicular to the chain was optimized
employing the B3LYP exchange-correlation
functional <cit.>, since this functional was shown 
in previous studies <cit.> to reasonably reproduce the bond length alternation of tPA 
observed in experiment <cit.>. 
For the geometry optimization a tight 
tier-2 basis set and a 1× 1× 20 𝐤-grid was used. 
Due to the computational overhead associated with additional vacuum
when computed using a PW basis, the vacuum around the tPA chain was reduced to at least 7 A
for calculations with VASP.
The CC and GW calculations involving FHI-aims were performed using the loc-NAO-VCC-nZ basis sets
developed by Zhang et al. <cit.>.
The basis set convergence for the band gap of tPA was
investigated on the EOM-CCSD level of theory, the results of which are shown in Table
<ref>.



As Table <ref> conveys, the loc-NAO-VCC-2Z basis set allows to converge the band gap to
 157 meV for a 1× 1× 6 𝐤-mesh and 141 meV for a 1× 1× 8 𝐤-mesh, essentially independently of the size of the BvK cell, which we consider to be sufficient for the present application as the finite-size error rather than the basis-set incompleteness error is at the center of this study. Hence, all calculations
involving FHI-aims will be performed using the loc-NAO-VCC-2Z basis set. 

In order to approximate the contribution from the singularity of the Coulomb potential,
a variation of the spherical truncation approach pioneered by Spencer and Alavi <cit.> is employed
in FHI-aims. In this cut-Coulomb approximation the long-ranged part of the Coulomb potential 
is made to decay fast to 0 beyond a set radius r_cut  by multiplying the original 1/r
potential with a complementary error function <cit.>


    v^cut-Coulomb = 
    	1/2|𝐫|erfc[η(|𝐫| - r_cut)]


with η as the inverse decay width.































§ RESULTS

































To determine the validity of the convergence rate of the finite-size error
for IP- and EA-EOM-CCSD, the findings will be presented as follows:
First, the fundamental properties of the EOM-CC structure factor that have been
mentioned in part already in Section <ref>, will be elucidated
by looking at the structure factor of the 3D LiH primitive cell repeated
periodically in one direction. Following that qualitative discussion,
the previously derived convergence rates of the EOM-CC band gap finite-size error will
be applied to a one-dimensional system, a chain of trans-Polyacetylene. Finally, the 
analytical extrapolation expression will
be applied to the finite-size convergence of
the G_0W_0 band gap, testing its validity outside of EOM-CC theory.



 §.§ The EOM-CC structure factor of LiH in one dimension



In order to relate the herein newly derived EOM-CC structure factor to the
ground-state analogue, the 
transition structure factor <cit.>, we
want to start by computing these two quantities for the LiH chain 
using supercells, which consist of increasingly many unit cells in a single direction.
Since the purpose of this section is 
a direct comparison of the correlation structure factors of the two methods, only the Γ-point
was sampled in reciprocal space.
As noted in Section <ref>, the derived behavior
of S^IP/EA(q) for q→ 0 (and equally so for the
transition structure factor) should hold independently of the
system's dimension.

We shall start by reviewing the properties of the ground-state CCSD transition
structure factor, which is shown in Figure <ref> for
the LiH chain for super cell sizes of  1× 1× 10,
1× 1× 50 and 1× 1× 100. Subsequently, a comparison to the EOM-CC
structure factor will be drawn.

In analogy to the EOM-CC structure factor, which upon contraction with the Coulomb potential
and integration over the reciprocal space yields 
the correlation contribution to the charged quasi-particle energies (see Equation <ref>), 
the ground-state correlation energy E^corr_0 can be obtained in the same way by means 
of the transition structure factor S(𝐪)


    E^corr_0 = ∫d𝐪  S(𝐪)v(𝐪).









Figure <ref> shows the medium- to long-range portion of
the transition structure factor in the direction of the chain, that is the z-direction. Even though,
S(𝐪) can be computed for longer 𝐪-vectors, corresponding to short-range
processes in real space, this is not relevant to the present discussion of the finite-size convergence.
The vertical lines in Figure <ref> show the magnitude of the smallest 𝐪-vector
that can be resolved in the respective super cell. That minimal 𝐪-vector is given by the minimal distance
of two 𝐤-points or (in a super cell formulation) by the smallest reciprocal lattice vector and corresponds
to the biggest real-space distance captured by the bvk cell. As can be observed in Figure <ref>,
the bigger the super cell size of LiH, the smaller does this minimal 𝐪-vector become and
the more long-range correlation information
does the transition structure factor contain.
The transition structure factor in Figure <ref> features
a minimum, at some material- and state-specific distance in reciprocal space.
This distance can be interpreted as a characteristic distance for the electronic correlation, as 
we expect the
contribution of the CCSD transition structure factor to the correlation energy to be maximal in the vicinity 
of its extremum.
In the specific case of the LiH chain, one finds that the minimum is at 1.4 Å^-1 for 
the ground state correlation energy, which corresponds to approximately a 1× 1× 2 supercell.
Another fundamental property of the transition structure factor of the ground-state is the fact
that for |𝐪|→0, S(𝐪)→0 with zero slope as is apparent in 
Figure <ref>.
This directly results from the sum-rule of the pair correlation function, which is the Fourier transform
of the transition structure factor.

The EOM-CC structure factor for charged excitations exhibits qualitative differences to the ground-state
transition structure factor. The structure factor for the first IP and EA in the EOM-CCSD framework for the
LiH chain is shown in Figure <ref>.



The EOM-CC structure factors feature a minimum as well, but it appears at 
a significantly smaller |𝐪| value of 0.5 Å^-1, 
which corresponds approximately to a 1× 1× 6 super cell. This reflects the
longer range of correlation effects in the EOM-CC case in comparison to the ground-state.
This comparison
clarifies that the problem of reaching the tdl for charged quasi-particle energies is substantially
more difficult than for the ground state case.
Another point of departure between ground state CC
and EOM-CC theory, is the value of the structure factor at 𝐪=0. 
While in the case of the ground state the S(𝐪=0) value is always 0,
in the case of EOM-CC we have previously derived that S(𝐪=0) is a finite negative value in our
formulation corresponding to the many-body character of the respective excitation.

As Figure <ref> demonstrates, all the properties discussed for the EOM-CC structure factor apply in a similar fashion to both the IP- and the EA-EOM-CCSD
case and the IP+EA case.
Let us stress once more that the IP+EA case corresponds to the electronic band gap in the presently used convention. It should be stressed at this point, that the EOM-CC structure factor that we have introduced
in this work only captures the correlation contribution to the IP and EA quasi-particle energies and to
the band gap, which are given by the many-body contributions of the EOM-CC equations as discussed
in Section <ref>. In order to obtain the full excitation energy, we compute and converge
the generally sizable single-body contributions (compare Equation <ref>) to the EOM-CCSD expectation value separately.

With respect to the previously derived asymptotic behavior of the EOM-CC structure factor for 𝐪→ 0, it becomes immediately apparent from the LiH EOM-CC structure factors, particularly for the 1× 1×
100 cell, that it does not seem to exhibit significant linear behavior at 𝐪=0.
Furthermore, S(𝐪=0) is relatively small.
Therefore, one can conclude for the anisotropic LiH super cell,
that both the constant |𝐪|^0 and the linear contribution |𝐪|^1 contribution of the EOM-CC structure
factor to Δ_FS^IP/EA is negligible.































































 §.§ The trans-Polyacetylene chain




  §.§.§ The EOM-CC structure factor of tPA


In order to verify the applicability
of the convergence rates, which were derived in the
long-wavelength limit (|𝐪|→ 0), we now investigate 
the IP- and EA-EOM-CCSD structure factors of an actual one-dimensional system, the tPA chain 
(note, that the previosuly studied LiH system, is a bulk material, however only extended in one direction). 
For that purpose, super cell sizes
of up to 1× 1× 32 were computed via IP- and EA-EOM-CCSD.
The IP- and EA-EOM-CCSD structure factors, and the structure factor corresponding to the correlation
energy contribution to the band gap, that is the sum of IP
and EA in the present convention, are illustrated in Figure <ref>.

Figure <ref> shows that
for the 1× 1× 32 super cell of the tPA chain, the minimum
of the EOM-CC structure factor is resolved and two data points to the
left of the minimum are obtained.
As discussed previously for the case of the LiH EOM-CC structure factor,
the analytical extrapolation expressions derived
for the asymptotic limit (|𝐪→ 0|) are strictly speaking only applicable for
system sizes for which the minimum of the structure factor
can be resolved. By applying a cubic spline interpolation
to the EOM-CC structure factors in Figure <ref>,
one finds that the minimum is located near  0.2 Å^-1, which corresponds
to approximately a 1× 1× 13 bvk cell. It is, however, important to note that the data point directly
to the right of the minimum in Figure <ref>, which corresponds to a 1× 1× 10 bvk cell,
is in close proximity to the minimum itself, so that the quadrature error resulting from extrapolation of
system sizes as small as 10 primitive cells in one direction can be assumed to be small.
Compared to the characteristic distance of 0.5 Å^-1 previously found for LiH,
this suggests that the length scale of the electronic correlation is more than doubled compared to LiH.
One possible explanation for that is the presence of long-ranged dispersion interactions, which are expected to be
more prominent in an organic compound like tPA than in an ionic compound with small ions like LiH.
Note, that in the case of the IP-EOM-CCSD structure factor in Figure <ref>, we observe
a second, shallow minimum at 0.8 Å^-1. Since this feature is located in the medium- to short
range region of the EOM-CC structure factor, this is most certainly an artefact of the comparably small PW
basis set, which had to be used to compute the 1× 1× 32 super cell of tPA. As, however, we are interested
in the long-range characteristics of the electronic correlation, this is not expected to have any effect on the
properties of the tPA EOM-CCSD structure factors discussed so far.





  §.§.§ Convergence of the correlation contribution


If the derived convergence rate in Equation
<ref> does indeed accurately model the
convergence of one-dimensional systems to the tdl, will be determined by
studying the convergence of the numerically determined band gap for trans-Polyacetylene (tPA).





In the pursuit to determine if some order of q of
S^IP/EA(q) in Equation <ref> is
dominant and if the derived convergence rate in Equation
<ref> can be reduced meaningfully to a
single leading-order contribution, a set of convergence rate models
are applied to the finite-size convergence of tPA in Figure
<ref>.  In particular, the full derived
model in Equation <ref> is fitted to the
calculated band gaps as a function of N_k along the direction of the
tPA chain. This model is compared to three other convergence rates,
which are the leading-order contributions originating from the q^0,
q^1 and q^2 contribution to the long-wavelength limit of the EOM-CC
structure factor, which are AN_k^-1ln(N_k^-1),
AN_k^-2ln(N_k^-1) and AN_k^-3ln(N_k^-1),
respectively, where A is determined by fitting to the numerical data. 
To prioritize the accurate description of the finite-size
convergence in the long-range limit, a modified least-square
cost function was used for fitting, where a data point corresponding
to a system size of N_k was weighted with
a factor of N_k^2.
In agreement with the computed
EOM-CC structure factor of tPA shown in Figure <ref>,
it was found that the derived convergence rates only describe
the numerical data in Figure <ref>
qualitatively well for system sizes that exceed some minimal number
of 𝐤-points N_k, so that the models in 
Figure <ref> are fitted to the
EOM-CCSD band gaps for N_k≥ 8. To remain consistent, all
following fits to FHI-aims band gaps, that is 
Figure <ref>,
<ref> and
<ref>, are shown for 
N_k≥ 4, but only fitted to data points with N_k≥ 8.

Undeniably, the q^0-order contribution to Equation
<ref> on its own, given by A N_k^-1lnN_k^-1,
fails entirely to model
the computed data in Figure <ref>,
underestimating the band gap in the tdl by almost  500 meV
relative to the full model. This lends credence to the aforementioned
estimation that for quasi-particle excitations with a minor many-body
character, these q^0 contributions would only play a minor role in
the description of the convergence to the tdl and is also consistent
with the small contribution of the |𝐪|^0 term of the
EOM-CC structure factor to the correlation energy in the
anisotropic LiH cell in Section <ref>. As a matter of fact,
the many-body contribution of the IP and the EA in tPA were found
to be both about 5%. Moreover, we stress that even if this contribution was
large, the employed Coulomb singularity treatment in FHI-aims would capture
this contribution to the finite-size error.
Similarly, the model corresponding to a EOM-CC structure factor linear
in 𝐪, given by AN_k^-2ln(N_k^-1),
converges visibly slower than the calculated data points
leading to an underestimate of the band gap of over 130 meV
compared to the full model.
This, in combination with the observation that the LiH EOM-CC structure factor
does not seem to exhibit any linear behavior for q→ 0, hints at the
possibility that in some systems – or possibly in general – there is no
or only a negligible linear contribution to S^IP/EA(𝐪) for q→ 0.
Instead, at least in the case of tPA, the q^2 contribution to the EOM-CC
structure factor appears to dictate the tdl convergence of the band gap.
Even though the AN_k^-3ln(N^-1) model notably deviates from the calculated
data for smaller N_k, it matches perfectly for N_k=8 to N_k=24,
lying virtually on top of the full model and yielding a band gap that
is underestimated by 43 meV.

This range of applicability is
slightly larger than the prior investigation of the EOM-CC structure factors
for tPA would suggest, were a minimum of 10 primitive cells or 𝐤-points
in one direction were found to be necessary. One likely explanation for that
minor disagreement is the utilization of different Coulomb potential approximations
in both codes as detailed in Section <ref>.

To conclusively verify the precision and accuracy of both the EOM-CCSD band gap data obtained via FHI-aims 
and the herein proposed extrapolation approach, 
the EOM-CCSD band gap calculations for supercells of size up to 1× 1× 16 were repeated using
one of the structures investigated by Windom et al.<cit.>. In that study, the fundamental
band gap of different tPA geometries was investigated on the
EOM-CCSD level of theory employing increasingly long
tPA oligomers with a cc-pVTZ basis set. 
For the comparison, the central C2H2 unit of one of these structures, 
in the original paper denoted by tPA3,
was extracted and treated under periodic boundary conditions
in the same manner as the B3LYP-optimized structure studied
so far.
To ensure a fair comparison, the  calculations
were performed employing the loc-NAO-VCC-3Z basis, the results
of which are shown in Figure <ref>.



Using the previously identified AN_k^-3ln(N_k^-1)
leading-order model, an EOM-CCSD band gap of 4.914 eV
in the bulk-limit is found.
This is in reasonable agreement
with the estimated tdl value of 5.07 eV of the original study. The
remaining deviation of roughly 160 meV is likely the
result of comparably small oligomer sizes, namely 6 to 9 C2H2 units 
that were used for the extrapolation. Also, the fact that
neither of the results are
converged with respect to the basis set (see Table <ref>) makes an exact
comparison more difficult.










  §.§.§ Convergence of the mean-field contributions


What has been neglected so far, however, is the potential influence of the single-body contributions
from the underlying hf calculations on the overall
convergence of the band gap. Most importantly,
the convergence of the hf exchange with
respect to system size needs to be considered.
So far, the hf calculations have been performed with the same system sizes as were subsequently
used for the CC calculations, so that the total finite-size error of the final EOM-CCSD band gap contained
contributions from both the mean-field calculation and the correlated CC calculation.
To allow for a more systematic study, which only includes long-range correlation effects, the EOM-CCSD band gap results will be recomputed using converged hf single-particle states and eigenenergies from supercells of size 1× 1× 48 or more. Via this down-sampling approach, one can independently analyze the
convergence of the correlation energy of the EOM-CCSD band gaps. The finite-size convergence including the leading-order fits as they have been
shown in Figure <ref> is presented in Figure <ref>
for the down-sampled data.



The fitted leading-order models in Figure <ref> suggest
that after the single-body contributions, particularly the hf exchange, has been accounted for, still the AN_k^-3ln(N_k^-1) model which results from the contribution to the structure
factor quadratic in q, describes the convergence
to the thermodynamic limit best and virtually lies on top of the data points for N_k > 10.



 §.§ Comparison to G_0W_0

To test the validity of the derived convergence rate for one-dimensional systems
outside of CC theory, the leading-order contribution determined in the previous
section is applied to the GW approximation, as well.
For that purpose, for the same system sizes that were previously computed via
IP- and EA-EOM-CCSD, the band gaps were obtained using the G_0W_0 method
with a hf starting point (G_0W_0@HF) to allow for a direct comparison to
the EOM-CC results.



The performance of the previously extracted AN_k^-3ln(N_k^-1)
convergence rate applied to both EOM-CCSD and G_0W_0@HF is shown in
Figure <ref>. There, one observes that
the convergence behavior of the GW method with respect to the number
of 𝐤-points N_k almost mirrors the one of the EOM-CCSD
method. In the same way the leading-order term originating from the
q^2 contribution to the EOM-CC structure factor models the band gap
convergence of both theories accurately over the entire range of
N_k=8 to N_k=24.

In an analogous step as in Section <ref>, the comparison
between EOM-CCSD and G_0W_0 band gaps was repeated after converging
the underlying hf orbitals and eigenenergies via down-sampling. As previously, for
that purpose the hf calculation was performed on a 𝐤-grid
of size 1× 1× 48 or more and the resulting hf orbitals and eigenergies
were subsequently used to compute the IP/EA-EOM-CCSD and G_0W_0
band gaps on 𝐤-meshes between 1× 1× 8 and 1× 1× 24.
The resulting convergence to the bulk-limit for the two methods
is shown in Figure <ref>.




Figure <ref> shows
that after the convergence of the single-body contributions of
the underlying hf method has been ensured, the remaining
contributions to both the EOM-CCSD and G_0W_0 quasi-particle energies
converge with the previously identified leading-order model
AN_k^-3ln(N_k^-1). While the down-sampling reduces
the individual band gaps of both methods most notably for small values
of N_k, the convergence behavior and the band gap in the
tdl remain unaffected, while the extrapolated G_0W_0@HF
band gap in the tdl changes by less than 6 meV

Figures
<ref> and
<ref>
also show the differences between 
G_0W_0 and EOM-CCSD band gaps retrieved as a function of N_k without and with applying
the down-sampling technique, respectively.
We note that the convergence of the difference can also be well
approximated using an AN_k^-3ln(N_k^-1) model.
The relatively small fitting parameter indicates
that both methods capture contributions to the
correlation energy with a similar magnitude and the same leading-order behaviour.
In other words, although the G_0W_0 and EOM-CCSD band gaps exhibit a
finite size error that is similar in magnitude and converges with the same analytical
behavior to the tdl, there exist
different leading-order diagrammatic contributions of both methods.
In practice one can still benefit from the similarity of both methods
by correcting the finite size error of EOM-CCSD band gaps using
the computationally significantly cheaper G_0W_0 method.
As noted before, the long-range behavior of
the EOM-CC structure factor is independent of the dimension of the
system. Therefore, one can infer that a GW-aided finite-size
correction technique for EOM-CC theory can in principle be performed
for one-, two- and three-dimensional systems alike.

An important question remaining is: What diagrammatic contributions
account for the differences between G_0W_0 and EOM-CCSD in the long-wavelength limit?
Although, the present work does not provide an explicit answer,
we note that previous studies performed detailed investigations of the 
relationship between the GW approximation and the
EOM-CC framework.
In the first study of this sort, performed by 
Lange and Berkelbach<cit.>, it was found that
the G_0W_0@HF approximation features more higher-order ring
diagrams than the Green's function of IP- and EA-EOM-CCSD.
These diagrams are known to be particularly relevant
for the description of electronic correlation in the long-wavelength limit of the ground state
for metallic systems <cit.>.
Since then, different flavors of the GW approximation have been
reformulated in a CC-like fashion. Bintrim and Berkelbach <cit.>
presented the EOM-CC-like working equations for the G_0W_0 method
in the Tamm-Dancoff approximation (G_0W_0-TDA) resulting in a frequency-independent
method.
More recently, the exact equivalence between the G_0W_0@HF method and the unitary
IP/EA-EOM-CCSD method using the quasi-boson approximation has been derived <cit.>.
We also note that different Green's function methods have been formulated based 
on the CC formalism <cit.>,
offering an alternative avenue to compare CC and GW methods.
The above relationships shall be further analysed in future work to develop computationally efficient and
accurate finite-size corrections for EOM-CCSD band gaps.



§ CONCLUSION

We have investigated the convergence of the EOM-CC band gap to the tdl
by a formal analysis of the correlation structure factor of the IP- and
EA-EOM-CCSD method S^IP/EA(𝐪) in the long-wavelength limit
(|𝐪|→ 0). As a result, we derived the convergence rate for
the one-, two- and three-dimensional case.  In order to verify the validity of
that approach, we focused on the one-dimensional case to be able to
compare to numerical results, using a chain of LiH unit cells and a
chain of trans-Polyacetylene as examples. Visualizing the
EOM-CC structure factor and modeling the finite-size
convergence using the derived convergence rate for one-dimensional systems suggests
that the band gap converges to the tdl with AN_k^-3ln(N_k^-1)
in leading order for N_k →∞. 
In analogy to the one-dimensional case, we expect a AN_k^-3/2 
behavior for two-dimensional systems and a AN_k^-1 behavior 
for three-dimensional ones.
Finally, we verified that our findings extend
beyond EOM-CC theory and apply to the G_0W_0 method as well.
We find that band gaps converge with the same rate in both
theories, 
providing a formal justification to use the GW approach for the extrapolation
of single-k-point EOM-CC results to the tdl.





  E.M. is thankful to Min-Ye Zhang and Sebastian Kokott for valuable discussions. 
  This project was supported by TEC1p [the European Research Council (ERC) Horizon 2020 research 
  and innovation program, Grant Agreement No.740233].  
  
  
  








































*



