Direct searches for general dark matter-electron interactions with graphene detectors: Part I. Electronic structure calculations

We develop a formalism to describe electron ejections from graphene-like targets by dark matter (DM) scattering for general forms of scalar and spin 1/2 DM-electron interactions and compare their applicability and accuracy within the density functional theory (DFT) and tight binding (TB) approaches. This formalism allows for accurate prediction of the daily modulation signal expected from DM in upcoming direct detection experiments employing graphene sheets as the target material. A key result is that the physics of the graphene sheet and that of the DM and the ejected electron factorise, allowing for the rate of ejections from all forms of DM to be obtained with a single graphene response function. We perform a comparison between the TB and DFT approaches to modeling the initial state electronic wavefunction within this framework, with DFT emerging as the more self-consistent and reliable choice due to the challenges in the embedding of an appropriate atomic contribution into the TB approach.


I. INTRODUCTION
Dark Matter (DM) plays a key role in simultaneously explaining otherwise anomalous physical phenomena that occur on extremely different astronomical length scales [1].For example, it provides the initial density fluctuations that trigger the formation of cosmic structures and generate the anisotropy pattern observed in the cosmic microwave background temperature and polarisation maps [2].It bends the light emitted by distant astrophysical sources, giving rise to spectacular gravitational lensing events [3] and, furthermore, it provides the mass required to support the flat rotation curves of spiral galaxies [4].
Despite these remarkable observations, we still do not know what DM is made of.The leading hypothesis in astroparticle physics is that DM is made of unidentified, yet-to-be-discovered particles [1].While this simple assumption can collectively explain all phenomena listed above, the hypothetical particles forming our universe's DM component have so far escaped detection.This dichotomy between solid gravitational evidence and lack of microscopic description makes the search for the "DM particle" a top priority.
A prominent class of experiments searching for the DM particle relies on the direct detection technique [5,6].This technique seeks for rare interactions between DM particles from the Milky Way and detector materials located deep underground in low background environments.As far as the DM-material interaction is concerned, DM direct detection experiments have until recently focused on the search for nuclear recoil events induced by the scattering of Weakly Interacting Massive Particles (WIMPs) in crystals, or liquid noble gases [7].Consequently, direct detection experiments have so far only probed DM particles of mass above about 1 GeV, as lighter particles would not be able to cause an observable nuclear recoil.However, the lack of detection of WIMPs has recently motivated the exploration of alternative experimental approaches that are better suited to probe DM particles of sub-GeV mass [8].It is in this context that DM direct detection experiments sensitive to DM-induced electronic transitions or electron ejections in materials play a central role.
Materials that have been proposed to search for sub-GeV DM particles via DM-electron interactions include liquid argon [9] and xenon [10][11][12], semiconductor crystals [13][14][15][16][17][18][19][20][21][22], 3D Dirac materials [23,24], graphene [25,26] and carbon nanotubes [27][28][29][30], to name a few.In this context, anisotropic media, particularly materials with anisotropic Fermi velocities such as graphene and carbon nanotubes, are interesting, as the associated rate of DM-induced electron ejections exhibits an enhanced daily modulation.This enhancement is caused by the structural anisotropy of the target material in combination with its relative orientation to the DM wind.Given that such a daily modulation is not present in typical experimental backgrounds, it would thus be a smoking gun for DM signal.A proposed experiment to search for DM-induced electron ejections from graphene sheets or arrays of carbon nanotubes, which is currently in the conceptual design stage, is the Princeton Tritium Observatory for Light, Early-Universe, Massive-Neutrino Yield, or PTOLEMY [31][32][33].
PTOLEMY's experimental design employs a large area surface-deposition tritium target coupled to a graphene substrate to detect the cosmic neutrino background via the observation of single electrons produced in the neu-trino absorption by tritium atoms [31].The coupling between tritium and graphene reduces the energy dispersion of the final state electrons by about an order of magnitude compared to the case of molecular tritium [32].A small electron energy dispersion allows for better discrimination between electrons produced by neutrino absorption and electrons populating the tail of the tritium β-decay spectrum [32].For this experimental setup to work, it is crucial to experimentally validate the use of graphene as a substrate by accurately measuring the electron-graphene interaction properties.In this intermediate stage of the PTOLEMY experimental program, when tritium target and graphene substrate are still decoupled, PTOLEMY can also operate as a directional MeV-scale DM detector.Specifically, two experimental configurations have been proposed.In the first one, a sample of stacked graphene sheets is considered (PTOLEMY-G 3 ) [25].Once an electron is ejected from one of the graphene sheets, it drifts in an external electric field until it reaches a calorimeter at the edge of the detector volume.This configuration allows for a full reconstruction of the final state electron kinematics.In a second experimental configuration (PTOLEMY-CNT) [27][28][29][30]33], an array of single-or multi-wall metallic carbon nanotubes is positioned in vacuum.When an electron is ejected from one of the nanotubes, it is driven by an electric field to the detection region, and recorded by a single electron sensor.The idea of adopting graphene sheets and carbon nanotubes as targets for directional, light DM detection has recently been further developed by the "Graphene-FET" and "dark-PMT" projects, respectively [34].
The possibility of using graphene or carbon nanotubes as directional detectors sensitive to DM-induced electron ejections motivates an accurate and comprehensive modelling of DM scattering by electrons bound in this class of anisotropic materials.In contrast, the rate of DMinduced electron ejections from graphene sheets [25] and from carbon nanotubes [27][28][29][30] has so far been computed assuming that the amplitude for DM-electron interactions depends on the momentum transfer only, and that the DM-electron interaction is spin-independent.This is a rather restrictive assumption, which can easily be violated, e.g. in models where DM has a non negligible magnetic or anapole moment [35].Furthermore, current estimates rely on the tight-binding approximation and have not been validated against first-principles calculations.
The purpose of this work is to extend and improve the formalism currently used to model the scattering of DM particles by electrons bound to graphene sheets.First, we extend the current formalism to virtually arbitrary DM-electron interactions by using the non-relativistic effective theory framework we developed in [35,36], and recently applied by the XENON group in an analysis of the electron recoil data reported in [37].Second, we improve the existing formalism by performing state-of-the art density functional theory (DFT) calculations in order to accurately model the electronic properties of graphene.
We expect that the formalism and findings we present here will be useful in the design of the PTOLEMY detector, as well as for the development of the Graphene-FET and dark-PMT projects.However, the relevance of our formalism goes beyond its application to these experimental concepts, as it can also be straightforwardly used to study the ejection of electrons in other experimental settings, where the final state is a free electron that can be described by a plane wave.We leave this exploration for future work.
This paper is the first of a two part series studying DM-electron scatterings in graphene targets.In this paper (Paper I), we lay the theoretical foundations.In the companion Paper II, we will focus on more explicit experimental setups and sensitivity studies [38].In addition, the software tools Darphene and QEdark-EFT developed for TB and DFT calculations respectively are publicly available [39,40].
This article is organized as follows.In Sec.II we introduce our general formalism for modeling the ejection of electrons by the scattering of DM particles in twoand three-dimensional periodic systems.In Sec.III, we describe the detailed electronic structure calculations we performed for graphene, both within the tight-bindingapproximation and within DFT.We apply these results to study the daily modulation of the DM-induced electron ejection rate for a hypothetical graphene detector in Sec.IV and conclude in Sec.V. We complement this work with appendices where we provide analytic formulae that are useful for evaluating our general electron ejection rates.

II. RATE OF ELECTRON EJECTION CAUSED BY GENERAL DARK MATTER-ELECTRON INTERACTIONS
In this section, we derive an expression for the rate of electron ejection caused by general DM-electron interactions in periodic systems.In Sec.III, we will perform the detailed electronic structure calculations that will enable us in Sec.IV to apply this general expression to the specific, experimentally relevant case of graphene.

A. General formalism
We are interested in processes in which a DM particle χ of mass m χ , initial velocity in the detector rest frame v, and momentum p = m χ v is scattered by an electron in initial state |e 1 .During the interaction, the DM particle transfers momentum q = p − p to the electron, where p is the final DM momentum, and causes an electronic transition from |e 1 to the final state |e 2 .In the notation of [35,41], the rate R 1→2 for this transition is where m e is the electron mass, n χ = ρ χ /m χ is the local DM number density, ρ χ = 0.4 GeV cm −3 is the local DM mass density, and f χ (v) is the local DM velocity distribution boosted to the detector rest frame.For f χ (v), we assume a truncated Maxwell-Boltzmann distribution, as in the so-called Standard Halo Model (SHM) [42].Specifically, and we take v 0 = |v 0 | = 238 km s −1 [42] for the local standard of rest speed, and v esc = 544 km s −1 [42] for the galactic escape speed.Following [24], we express the Earth's velocity with respect to the galactic centre, v e , in a coordinate system with z-axis in the v 0 + v direction, v the Sun's peculiar velocity and v e = |v 0 + v | 250.5 km s −1 [42], sin α e sin β sin α e cos α e (cos β − 1) cos 2 α e + sin 2 α e cos β

  ,
where α e = 42 • , β = 2π t/day, and t is the time variable.Finally, we also introduced the normalization constant, The total initial (final) energy E i (E f ) in Eq. ( 1) is the sum of the DM and electronic energies, where E 1 (E 2 ) is the energy eigenvalue of the electronic state |e 1 (|e 2 ).We denote the corresponding wave functions by ψ 1 and ψ 2 , and their associated Fourier transforms by ψ 1 and ψ 2 , respectively.These electron wave functions enter the electron transition amplitude M 1→2 , defined as in Eq. ( 14) of [35] by the integral, where M( , p, q) is the free electron scattering amplitude, and the initial state electron momentum.Here, we use momentum conservation to eliminate explicit dependence on the final state electron momentum from M. Furthermore, since the scattering of Milky Way DM particles by free electrons is expected to be nonrelativistic, we use the Galilean invariance of M to write M( , p, q) = M(q, v ⊥ el ), where v ⊥ el = v−q/(2µ χe )− /m e and µ χe is the DM-electron reduced mass.Finally, we expand M at linear order in /m e , and write it as follows [35] This expansion allows us to express the transition amplitude as where we introduce the scalar and vectorial overlap integrals, In order to evaluate the expressions above, we need to specify the initial and final electron wave functions.We begin by specifying the final-state electron wave function for the case in which the electron is ejected by the DM particle.In this case, the state |e 2 asymptotically approaches a free particle of momentum k .Consequently, the wave function ψ 2 (x) can be approximated by the plane wave which is normalised to one over a finite volume V .For electrons initially bound in graphene, this plane wave approximation has been validated by comparing results from angular-resolved photoemission spectroscopy (ARPES) measurements with simulated photoemission intensity maps, for which excellent agreement was found [43].Within this plane-wave assumption, we can express the scalar and vectorial overlap integrals in Eq. ( 8) and Eq. ( 9) in terms of the Fourier transform of the initial state electron wave function Also, using a plane wave as a final state, we find that the square of the transition amplitude in Eq. ( 1) can be written as where we introduced the free particle response function R free (k , q, v), for which we give a general expression in Appendix A. In order to understand the physical meaning of R free (k , q, v), it is instructive to take the limit of a free initial state electron in Eq. ( 13), and hence replace ψ 1 (x) with a plane wave of linear momentum .In this limit, one finds where all information, besides momentum conservation, is contained in R free .This shows that the second factor in Eq. ( 13) contributes non trivially to the squared transition amplitude only when the initial state electron is bound within a material.In this latter case, it encodes all relevant material properties via the Fourier transform of the initial state electron wave function.As we will see in the next section, this factorization allows us to express the rate of DM-induced electron ejection from materials in terms of a single material response function.This is in contrast with our previous findings for the cases of atomic ionizations [35] and excitations in crystals [41] where up to five material response functions were required to evaluate the rate of DM-induced electronic transitions between filled valence and empty conduction bands.One should also note that the results reported in [35,41] neglect the directional information of the event rate and assume a simplified treatment of the velocity integral in the transition-rate formula.By performing this integral exactly, as we do here via Monte Carlo integration (see Sec. IV), up to five scalar and two vectorial material response functions are in general expected to contribute to the DM-induced electronic transition rate [41].

B. Effective theory expansion of the scattering amplitude
In order to evaluate our general electron ejection formulae, Eqs.(1) and (13), we need to specify the coefficients, M(q, v ⊥ el ) =0 and ∇ M(q, v ⊥ el ) =0 in the nonrelativistic expansion of the scattering amplitude M in Eq. (7).From these coefficients, one can in turn obtain an explicit expression for the free-particle response function R free , as shown in App. A. In this work, we calculate these coefficients using effective theory methods.Specifically, we extract them from the non-relativistic effective theory of spin 0 and spin 1/2 DM-electron interactions [35], within which the scattering amplitude can be written as Here c i is the dimensionless effective coupling corresponding to the i-th operator, O i , in Tab.I, angle brackets denote an expectation value between DM-electron spin states, and F DM,i (q) is the DM form factor that encapsulates the q-dependence of the amplitude not captured by the operator O i itself [13] 1 .Possible forms of the DM form factor include 1 for short-range interactions, for long-range interactions, where we introduced an arbitrary reference momentum transfer q ref .In the context of sub-GeV DM searches, q ref is usually set to αm e , where α is the fine-structure constant, since this is the typical momentum of an electron in the outer atomic orbitals.Eq. (15) gives the most general expression for the nonrelativistic amplitude for DM-electron scattering that is compatible with momentum conservation and Galilean invariance.The formalism we develop in this work, and in particular the free-particle response function used in the numerical calculations, relies on the expansion in Eq. (15).

C. Benchmark particle physics models
With Eq. ( 15), we have a general parametrization of the non-relativistic scattering amplitude that virtually any fundamental DM particle model can be mapped onto.In the numerical applications presented in Sec.IV, we focus on four benchmark models, each of them corresponding to a different linear combination of operators 1 If the non-relativistic amplitude M(q) contains a given operator O i within two terms with distinct q-dependencies, i.e. two different DM form factors, one can still use our formalism by replacing c i F DM,i (q) with the sum c i F DM,i (q) + c DM,i (q) for that particular operator.
Interaction operators defining the non-relativistic effective theory of spin 0 and 1/2 DM-electron interactions [35].Se (Sχ) is the electron (DM) spin, v ⊥ el = v− /me − q/(2µχe), where µχe is the DM-electron reduced mass, v ⊥ el is the transverse relative velocity and 1χe is the identity in the DM-electron spin space.In the case of elastic scattering, v ⊥ el • q = 0, which explains the notation.
in Tab.I.These models -briefly reviewed below -provide interesting examples of DM-electron interactions, and demonstrate both the generality of the effective theory expansion in Eq. ( 15), as well as how the mapping from fundamental to effective coupling constants works in practice.

Dark photon model
Our first benchmark model has guided the direct search for sub-GeV DM particles over the past few years, and is referred to as the dark photon model.In this framework, the Standard Model (SM) Lagrangian is extended by a new U (1) gauge group with gauge coupling g D , and by a massive dark photon A [13,[44][45][46][47][48].The DMordinary matter interaction portal opens via a kinetic mixing between ordinary and dark photons in the interaction Lagrangian, i.e. εF µν F µν , where F µν (F µν ) is the field strength tensor of the ordinary photon (massive dark photon).The Lagrangian of the dark sector in this model is given by with the covariant derivative defined as where g D is the gauge coupling corresponding to the dark U (1) gauge group.
In our general framework, the DM-electron scattering amplitude in the dark photon model can be mapped on to the operator O 1 if one relates the coupling constants g D and ε to the effective coupling c 1 as follows, with

Electric dipole interactions
As another less trivial example that illustrates the wide applicability of our framework, we next consider the case of electric dipole DM-electron interactions induced by the interaction Lagrangian In [35], we showed that the scattering amplitude of this model can be mapped to the operator O 11 via

Magnetic dipole interactions
Similarly, one can assume an interaction portal via magnetic dipole interactions between DM and electrons by the following interaction term in the Lagrangian, As shown in [35], the corresponding scattering amplitude can be identified with a linear combination of four of the operators in Tab.I, with non-zero effective couplings given by As one can see, in this case the amplitude M is a linear combination of "short-range" and "long-range" contributions.

Anapole interactions
Finally, we also study anapole interactions, defined by the interaction Lagrangian Just as before, we compare the scattering amplitude with the effective operators and find a correspondence to O 8 and O 9 with the effective couplings with F DM,9 (q) = 1 .(25b)

D. Application to periodic systems
We continue the development of our formalism by specifying the wave function for the initial state electrons, which we assume to be bound within a crystal consisting of periodically repeating atoms.In this case, the electron eigenfunction is characterized by an energy band, i, and a lattice momentum, k, and has the form of a Bloch state, in which with a being an arbitrary lattice vector.We refer to App.C 1 for further details on Bloch's theorem and Bloch states.In this work, we consider two different bases to form the Bloch states of the initial state electrons.For our tight-binding approximation calculations (Sec.II D 1) we build the wavefunctions from linear combinations of atomic orbitals, and for our DFT calculations (Sec.II D 2) we use linear combinations of plane waves.
Assuming a Bloch state of crystal momentum k for the initial electron wave function and a plane wave of linear momentum k for the final state electron, we can combine Eq. ( 13) with Eq. ( 1) to calculate the transition rate R ik→k , that is the rate of electron ejections by DM scattering when the initial electron is in energy band i with crystal momentum k and the outgoing electron has linear momentum k .The energy difference in Eq. ( 1) now reads where E i (k) is the energy of an electron in energy band i with wavevector k, which is negative for bound electrons.Φ is the (positive) work function, which corresponds to the energy difference between the highest occupied electronic state and the zero-energy unbound free-electron plane-wave state; in this work we take the measured value for graphene of Φ = 4.3 eV [25,43].
Adding the contributions from all occupied initial states yields where the factor of 2 accounts for the double occupation of each electronic state due to spin degeneracy.We then sum the contributions from all final plane-wave states to obtain the total rate of electron ejections by DM scattering, Finally, we express the total ejection rate, R, as follows where we define the electron's energy change as ∆E e ≡ k 2 2me +Φ−E e , and introduce the material-specific response function As implied by the above δ-function, E e takes the value of the initial state electron energy where E e = 0 eV for the highest occupied state and E e < 0 eV for other bound states.Eq. (31) shows that both free-electron physics and material properties contribute to the total rate of electron ejections by DM scattering in a factorisable way, and that this factorisation involves a single material response function.The response function is normalized to where N bands is the number of occupied initial-state electronic bands and we use the normalization of ψ ik ( ) (see Eq. (D3)) and In the next two subsections we express the response function W by writing ψ ik as a linear combination of atomic orbitals and plane waves.Atomic orbitals are often employed within the tight-binding approximation (see Sec. III A), whereas plane waves are a standard basis in DFT electronic structure calculations (see Sec. III B).

Atomic orbital basis
Let us now derive a compact expression for the response function W by expanding the initial state electron wave function, ψ ik (x), in a basis of atomic orbitals, where j runs over the n atomic orbitals present in each unit cell, and Φ jk (x) are the Bloch states corresponding to each atomic orbital, Here ϕ j is an atomic wave function on an atom at position R jr , N cell is the number of unit cells and N k is a normalisation constant defined in App.C 1. Within the tight-binding approximation (introduced below in Sec.III A), the coefficients C ij in Eq. ( 34) are computed by solving the secular equation (Eq.C12) with band energies extracted from measurements, or calculated using a more sophisticated technique such as DFT.
To evaluate the material response function defined in Eq. ( 32), W ( , E e ), we need to calculate the square of the Fourier transform of ψ ik (x).Denoting by ϕ j the Fourier transform of ϕ j , for the Fourier transform of ψ ik we find The lattice sum can be evaluated by using Eq.(A.8) of [49].
where δ j is the location of the atom hosting the jth orbital in the unit cell that contains the origin of coordinates.For a given j, δ j = 0 if there exists a r ∈ {1, . . ., N cell } such that R jr = 0, and δ j = 0 otherwise.In Eq. ( 36), the sum runs over the reciprocal lattice vectors G.For each j and r, they satisfy G • (R jr − δ j ) = 2πm, where m ∈ Z. Another useful identity to evaluate the squared modulus of ψ ik is where 36) and ( 35), we obtain and by combining Eq. ( 38) with Eq. ( 32), we obtain the following expression for the response function W : Here, we performed the integral over the lattice momentum k, such that the δ function fixes it to k = − G * , where G * is the unique reciprocal lattice vector that ensures that k lies within the first Brillouin zone.The other terms of the sum over the vectors G do not contribute.
To evaluate Eq. ( 39), one needs to specify the wave functions ϕ j , which we will do in Sec.III A.

Plane wave basis
Let us now express the response function W by using plane waves to write the electron wave function ψ ik (x) as where 1 for all i and k.From Eq. ( 40), we find and where we used V = (2π) 3 δ (3) (0).This result can be directly inserted into Eq.( 32), which leads to the response function In Sec.III B 3, we extract the u i (k + G) coefficients and the band structure E i (k) from state-of-the-art DFT calculations.

III. ELECTRONIC STRUCTURE CALCULATIONS FOR ELECTRON EJECTIONS IN GRAPHENE DETECTORS
The equations derived in the previous section refer to three-dimensional periodic systems.We are now interested in applying them to the specific case of graphene, which is a single-layer material that is periodic in two dimensions.The "dimensional reduction" can be performed straightforwardly by means of the replacement specified below, where A cell is the two-dimensional unit cell of graphene, while k in the right-(left-)hand-side is a two-dimensional (three-dimensional) lattice vector in the first Brillouin zone.
With a general formalism for electron ejections by DM scattering in graphene detectors in place, we can now focus on the evaluation of the predicted electron ejection rate.This crucially depends on the response function W , which is in turn a function of the initial state electron wave functions.As a result, numerical evaluation of the predicted electron ejection rate requires detailed electronic structure calculations for graphene.In the following, we perform such electronic structure calculations using two methods: the tight binding approximation, and DFT.From this analysis, DFT will emerge as our recommended framework for electronic structure calculations for DM-electron scattering in graphene-based DM detectors.

A. Tight binding
To obtain the graphene response function in the tight binding (TB) approximation, we need to evaluate Eq. (39).The missing ingredient at this point are the coefficients C ij (k) that yield the contribution of atomic orbital j in band i to the response function.We separate the π-and σ-electrons and write the response function as In the TB approximation, the coefficients are found as the eigenvectors of the generalized eigenvalue problem in Eq. (C11).For a detailed review of the TB approximation in general and for the specific case of graphene, we refer to App.C 1 and C 2 respectively.

π-electrons
In case of the π-electrons in graphene, this eigenvalue problem can be solved analytically, as described in App.C 2. Therein, the full wavefunction of the πelectrons are derived in position and momentum space, which can be found in Eq. (C24) and (C33).The eigenvalues E π (k) and eigenvectors C π , required for the response function, are given by Eqs.(C22) and (C23) respectively.Therefore, the π-electron contribution to the response function can be written out explicitly.It can be shown to simplify to Here, the phase ϕ k and the normalization constant N k are given by Eqs.(C23) and (C29) respectively.

σ-electrons
While the procedure for the σ-electrons is conceptionally identical, the fact that their wavefunctions involve combinations of three atomic orbitals at two atomic sites means that the generalized eigenvalue problem of Eq. (C11) involves 6 × 6 matrices, which can no longer be solved analytically.Instead, we rely on numerical procedures where we use the Eigen library [50].
The involved six-dimensional matrices H and S are listed in Eq.(C34) of App.C 2. Using the numerical procedures of the Eigen library, we obtain the eigenvalues or band energies E σi (k) as well as the eigenvectors C σi . 2 Finally this allows us to evaluate the σ-contribution to the response function, By adding up the four distributions given by Eqs. ( 46) and ( 47), we obtain the final TB estimate of the graphene response function.Finally, we point out that our TB treatment of the electrons in graphene differs from a previous treatment by Hochberg et al. [25], in particular with regard to the Bloch wavefunctions given by Eq. (34).A second crucial difference is our choice for the atomic wavefunctions, which we discuss next.We present a detailed comparison in App.E.

Atomic wavefunctions
In order to evaluate the TB estimates of the graphene response function given in Eqs. ( 46) and ( 47), we need to specify the atomic orbitals ϕ i (x) (or rather their Fourier transforms ϕ i ( )) for the electrons in carbon.
We expand the atomic wavefunction into a radial and directional component, , where r = |x| and Y m l (x) are spherical harmonics.Following [51], we describe the radial part of the atomic orbitals of carbon electrons as linear combinations of Slatertype orbitals (STOs). with We include a more detailed description of these wavefunctions, including the values of the different coefficients and the expressions for the Fourier transforms, in App.D 2. This choice differs from the previous approach by Hochberg et al. [25], who use re-scaled hydrogenic wavefunctions to approximate the electron wavefunctions in carbon atoms.We comment on this in App.E.

Limitations
As described in greater detail in App.C, we can reproduce the measured band structure of graphene by adjusting the overlap and transfer parameters of the TB approximation.However, in particular the overlap parameter, e.g. the parameter s in Eq. (C7), is given by the overlap integrals of atomic wavefunctions at neighboring atomic sites.For a given choice of atomic wavefunctions, it is therefore possible to compute s independently of the band structures.This gives rise to the issue of selfconsistency of this approach.
In [25], the authors use hydrogenic wavefunctions to approximate the atomic orbitals of carbon.As described in detail in App.D 1, they ensure consistency between the two independent values of the overlap parameters, by re-scaling the effective charge factor Z eff .However, the resulting wavefunctions do not resemble the atomic wavefunctions of carbon, as we describe in App.D, which is why we chose to use RHF wavefunctions.While it is possible to perform a similar re-scaling of the RHF wavefunctions to establish consistency with the overlap parameters listed in Tab.II, this generally modifies the wavefunctions to the extent that they no longer describe electrons in carbon atoms.Therefore, it seems to be a limitation of the TB approximation to reconcile the use of realistic atomic wavefunctions with the overlap parameters that reproduce the material's band structure in a fully self-consistent manner.This issue is not a characteristic of our specific treatment but rather a general feature of the TB approximation itself reflecting the phenomenological nature of this approximation.

B. Density Functional Theory
Having carefully described the features and limitations of the tight-binding approach, we now report on our DFT electronic structure calculations.We start with a brief review of the main assumptions underlying DFT in Sec.III B 1. In Sec.III B 2, we provide a general argument supporting DFT as a framework for electronic structure calculations in the case of DM-induced electron ejections by graphene targets.In Sec.III B 3, we describe the details of our specific DFT implementation in a modified version of the QEdark-EFT code.

Assumptions
Density functional theory (DFT) [52,53] (for a review see [54]) is a widely used method for calculating the ground-state electronic properties of materials.It allows for explicit treatment of the chemistry and crystal structure without the need for empirically determined input parameters, and provides well-tested and computationally affordable approximations for the manybody electron-electron interactions.In addition, it has been implemented in numerous publically or commercially available computer codes that are convenient to use.
The theory is based on the Hohenberg-Kohn theorem [52] which states that, for electrons in an external potential (in this case provided by the charged atomic nuclei), the total energy is a unique functional of the electron density, with the ground-state density being the one that minimizes the value of this functional.The electronic ground state charge density can therefore be obtained variationally.
In practice the charge density, n e (x), is written as a sum over the so-called Kohn-Sham wave functions, ψ i (x) [53] of a fictitious auxiliary system in which the electrons are non-interacting.This mapping enables convenient computational solution of the many-body Schrödinger equation at the expense of an inexact description of the quantum mechanical exchange and correlation terms.These have been obtained numerically using Quantum Monte Carlo for the homogeneous interacting electron gas [55], and a number of well-tested approximations exist that are appropriate for different material systems.An additional widely used and well-established approximation divides the electrons into valence electrons, which are treated explicitly within the DFT calculation, and low-energy core electrons, which are combined with the nuclei in the external potential.This pseudopotential approximation drastically reduces the computational expense and is chemically well founded, since the core electrons are not involved in chemical bonding and are only minimally modified in the solid.The choice of pseudopotential, and the numerical and implementational details of the DFT calculation for graphene performed here, are discussed further in Sec.III B 3.
We note that, while the Kohn-Sham wave functions and energies do not formally correspond to true singleelectron wave functions and energies (except for the highest occupied level, which provides the ionization energy), in practice their dispersion is usually in remarkable agreement with measured photoelectron spectra.The Kohn-Sham band structure is therefore often treated as a proxy for an effective single-particle band structure in a periodic solid.Since the Hohenberg-Kohn theorem describes only the ground-state electron density, however, this is particularly ill-founded for unoccupied conduction-band states.The methodology and results we present here, however, do not rely on any physical interpretation of the Kohn-Sham wavefunctions, since our final electron states are unbound plane-wave states, and, as we show next, our response function is primarily determined by the ground-state charge density rather than the ground-state wavefunctions.

Motivations
An important observation we can draw from our general electron ejection rate formula is that the graphene response function W is directly related to the ground state electron momentum density, ρ e , defined as the Fourier transform of the electron charge density.Indeed, an explicit second quantisation calculation allows us to write ρ e as where is the mean ground state occupation number density, while a † ik and a i k are second quantisation creation and annihilation operators associated with the ψ ik and ψ i k Bloch states, respectively 3 .In Eq. ( 49), the i = i diagonal term is directly proportional to the response function W .The i = i off-diagonal term describes band mixing effects arising from electron-electron interactions across different bands.While in general n ii (k) = 0 for i = i , off-diagonal contributions to ρ e are expected to be subleading in the case of graphene, where electron-electron interaction and correlation effects induce variations in the band energies E i (k) of at most a few %, see for example Fig. 1 in [56].
This is an important observation because it implies that our DFT predictions for W are only marginally affected by the lack of a clear interpretation for the individual Kohn-Sham states, which, in principle, is one of the limitations of a DFT approach.These states contribute to W mainly through one very specific combination, namely, the Fourier transform of the electron charge angle brackets denote an expectation value on the ground state, while ψ ik (r) is the Bloch wave function for the initial state electron.
density, which, by construction, is self-consistently computed in DFT.
We find this observation a solid argument in favour of DFT as a theoretical framework for computing the graphene response function W .This conclusion is also corroborated by the good agreement found between the measured and DFT-calculated graphene Compton profile [57], which is the longitudinal projection of the electron momentum density, and thus closely related to W .

Numerical implementation
For the numerical evaluation of the graphene responses in the DFT framework, the QuantumEspresso v.6.4.1 [58][59][60] code was used, and interfaced with QEdark-EFT [40], an extension to the previously established QEdark [16] package.Since QuantumEspresso uses a plane-wave basis with periodic boundary conditions, we simulated the graphene sheet as a system containing sheets separated by a large but finite distance, L z .
For the self-consistent calculations, we used the C.pbe-n-kjpaw psl.1.0.0.UPF pseudopotential provided with the QuantumEspresso package, which includes the 2s 2 and 2p 2 electrons in the valence configuration and treats the 1s 2 electrons as core.The minimal suggested energy cutoff for the plane-wave expansion for this pseudopotential is 40 Ry for the wave function and 326 Ry for the charge density.We chose much larger values-2000 Ry for the wave function cutoff and 16000 Ry for the charge density cutoff-since for the case of dark-matter induced excitations, we are interested in the high-momentum tails of the electronic wave functions that are usually unimportant for low-energy solid state physics applications.
We used the PBE exchange and correlation func-tional [61] with the experimentally measured lattice constant of a = 2.46 Å, and sampled reciprocal space with a 16 × 16 × 1 Monkhorst-Pack k-point grid, which is sufficient to capture the linearly dispersing Dirac cones at the Fermi level (see Fig. 5).Since carbon is a light atom relativistic effects are minimal and we did not including spin-orbit coupling.

Discretization
As a result of the widely spaced periodically repeating graphene sheets required by the periodic boundary conditions of our DFT code, Eq. ( 43) is discretised and the expression evaluated by QEdark-EFT becomes where δ E and δ are the bin size in energy and momentum respectively, and = k − q.Subscripts n, m, o and l denote the index of the corresponding momentum and energy bin, ω k are the weights of the reciprocal lattice k-points and, following the conventions of QuantumEspresso, k ω k = 2.The sum over i runs over 4 valence bands, and the sum over G is truncated by with the value of E cut = 27.2 keV.

C. Comparison of response functions
In Figs.1-4, we present different ways to visualize the graphene response function obtained with the TB and DFT methods.In these figures, we integrate the response function defined in Eq. ( 32) over all energies E e or an interval in E e to obtain a function of momentum only.
In Fig. 1, we depict the response function as a function of x and y with z ≈ 0. Both for TB and DFT, we find the characteristic hexagonal shape of the first Brillouin zone of graphene.In the case of TB, we also depict the individual contributions of the π and σ electrons.Note that the distortions of the DFT version originate in the finite grid sampling in reciprocal space and are washed out once we integrate over to obtain an observable.
In Fig. 2, we show the response function w.r.t. the angle between the initial state electron momentum and the orientation of the graphene sheet.The response function is integrated over all-electron energies as well as only over electrons within 5 eV of the Fermi level in order to understand the patterns observed in daily modulation plots for dark matter candidates of various masses.We can see that the DM particles carrying lower kinetic energies will be able to interact preferably with electrons that have their momenta pointing in θ ∼ 40 • whereas DM candidates with higher energy will be able to access electrons with momenta pointing in all directions.
In Fig. 3, we show the contribution of various electronbinding energies to the integrated response function for FIG. 5. Graphene band structure as calculated from TB and DFT.In the case of DFT, the Dirac cone crossing at Ksymmetry points was used to determine a precise value for the Fermi energy of the graphene sheet.We include the valence bands and the conduction band for the π-electrons.Further conduction bands are are not shown here, but can be found in Fig. 10 for TB.
both TB and DFT.Both approaches show a similar dependence of the momentum contribution integrated over all directions for the selected energy intervals.
In order to facilitate a quantitative comparison of the two approaches, we show the response function as a function of momentum ⊥ (perpendicular to the sheet) and (in-sheet momenta) in Fig. 4 using a log-scale on the y-axis.
Furthermore, the lower panel depicts the ratio of the W functions obtained by DFT and TB.We can see that for low momenta 15 keV, both response functions lie within a factor of 2 of each other.However, for larger momenta, the TB approximation predicts significantly larger values.
Overall, the Figs.1-4 demonstrate that TB and DFT predict the same qualitative features of the graphene response function.A more quantitative comparison reveals relative deviations between the two approaches that typically do not exceed a factor of 2. As mentioned above, the exception here is the case of large momenta, where we found larger deviations between the two methods.The response matrix at large momenta can be sensitive to the electron density close to the atomic nuclei, where the electron orbitals resemble the atomic orbitals the closest.This makes a qualitative comparison or assessment for large momenta difficult as these contributions are smoothened out in the DFT approach.FIG. 6. Calculated total time-averaged rate (left) and daily modulation curves (right) of DM-induced electron ejections from graphene obtained with DFT (blue) and TB (red).On the left, we show the rate for O1 contact type interaction (solid), O1 long range type interaction (dashed) and O3 contact type interaction (dash-dotted).The upper panel to the left gives the total rates as a function of DM mass, and the lower panel gives the ratio of the DFT to TB rates.On the right, we show the daily modulation curves for O1 contact type interaction, where the solid, dashed and dash-dotted curve is for a 5 MeV, 10 MeV and 100 MeV, respectively.

IV. CASE STUDY: DAILY MODULATION OF THE ELECTRON EJECTION RATE IN A GRAPHENE DETECTOR
After our study of the electronic structure calculations of graphene targets in Sec.III, we now present and compare the expected electron ejection rates that we obtain using both TB and DFT.We focus on a hypothetical experimental setting where the electron ejected by an incoming DM particle is recorded independently of the direction of ejection.We refer to a companion paper (from now onward, Paper II [38]) for detailed sensitivity studies of different settings for graphene-based DM detectors that are currently in a research and development stage.
Fig. 6 (left) shows a comparison of the time-averaged TB and DFT rates as a function of the DM particle mass for DM-electron interactions described by the operators O 1 (both contact and long range interactions) and O 3 (contact only) in Tab.I.
In the case of O 1 , we find that the electron ejection rates predicted by TB and DFT differ by less than a factor of 2 for most DM masses, and up to a factor of 3 at very low masses.Generally, we find that DFT predicts higher rates at low masses, whereas TB predicts higher rates at high masses.
In the case of O 3 contact type interactions, the quantitive comparison for low masses (m χ 20 MeV) is similar to O 1 .We find larger deviations for heavier masses, however, with the TB approach predicting an O( 10) larger rate at m χ = 100 MeV.For this particular operator, large momentum transfers are favoured, and these become kinematically more accessible for larger masses.
The difference in rate therefore originates in the graphene response function at large momentum , where TB predicts a higher response than DFT as seen in Fig. 4. The structure of O 3 therefore suppresses the contribution of the response function at low momentum, where we have better agreement between the two approaches.Here, the different treatment of the electron density close to the atomic nuclei plays an important role and obstructs a qualitative comparison.While this is the case of largest deviation between TB and DFT that we present, one should also note, that O 3 is an extreme but instructive case that does not arise from relativistic DM models at leading order [36].
Again comparing the DFT and TB computational frameworks, Fig. 6 (right) shows the expected daily modulation in the rate of DM-induced electron ejections for three values of the DM particle mass, and for the interaction operator O 1 .In all cases the expected electron ejection rate is divided by the corresponding timeaveraged rate, in order to cancel out the O(1) differences between the DFT and TB rates reported in the left panel of Fig. 6.From Fig. 6 (right), we conclude that irrespective of the chosen computational framework, the largest relative rate is found when the direction of the DM wind is perpendicular to the graphene sheet.This is due to the fact that the electrons are more spatially constrained in the out-of-plane direction, giving rise to higher momentum contribution to the electron wavefunction and thus increasing the total interaction rate (as discussed further in Paper II).The strong daily modulation we find for the rate of DM-induced electron ejections demonstrates that graphene is well suited for establishing a daily modulation signal characteristic of the directionality of the DM  wind.
In Fig. 7 we show the daily modulation pattern of electron ejections from graphene for various DM masses and interaction types, now focusing on DFT as a computational framework as the corresponding results from TB are qualitatively similar.The daily modulation pattern is similar for most of the DM masses and interactions with a maximum at around time=0h (when the DM wind is perpendicular to the graphene sheet) and a minimum at around time=12h (as for Fig. 6 (right)).However, for m χ = 2 MeV, the maximum is shifted to two peaks around time=4h and time=20h.This is a consequence of the 2 MeV DM particle only being able to eject electrons with E e close to 0, corresponding to the top two panels in Fig. 2. The location of the peaks around time=4h and time=20h is due to the DM wind aligning with the peak at θ ∼ 30 • and ∼ 4 keV.

V. SUMMARY AND CONCLUSION
In this paper, we have investigated two solid-statephysics approaches to modeling DM interactions with graphene-like targets, TB and DFT.Below, we summarise the main features of the two methods and the arguments that led us to identify DFT as the preferred framework for modeling the scattering of DM particles by electrons bound in graphene.
Both DFT and TB capture the main characteristics of the band structure of graphene, such as the Dirac cone at the K-symmetry point of the Brillouin zone and valence band energy distributions (Fig. 5), and predict a response function that reflects the symmetry of the reciprocal space underlying the graphene lattice (Fig. 1).
The two computational frameworks also predict a qualitatively similar daily modulation pattern in the total rate of electron ejections.However, the two methods generally predict electron ejection rates that differ by an O(1) factor with TB (DFT) predicting higher rates for high (low) DM masses.
The two approaches employ a different set of approximations that limit their predictivity and region of validity.In order to effectively perform DFT calculations, one usually chooses a radial cutoff to the pseudo-potential and smoothens the core-electron wavefunctions closer to the atomic nucleus.This approximates the electronic structure below that cutoff and suppresses some of the high-momentum contributions to the electronic wavefunction.The total event rate is therefore suppressed when higher energy excitations (of tens of eV) are considered [62].
On the other hand, DFT is able to provide a selfconsistent calculation of the ground-state electron density, which is the indirect quantity of interest for the DM interactions considered in this work.Indeed, an important result of this work is our demonstration that in cases where the outgoing electrons can be treated as a plane wave, the DM and electronic contributions to the DMinduced electron ejection rate factorize.The five crystal response functions identified previously in our work [41] then simplify into a single crystal response that is directly proportional to "diagonal part" of the Fourier transform of the ground-state electron density.
The TB approach has the advantage of easier implementability and computational affordability.In its usual low-energy applications, the detailed form of the atomic basis wavefunctions is not explicitly considered.However, since for the case of DM-electron scattering, we are interested in the explicit form of the electronic wavefunctions, one needs to embed the atomic wavefunctions into this framework.
This proves problematic since these atomic wavefunctions are required to satisfy the overlap integrals of TB that reproduce the experimentally observed band structure.The TB approach assumes that a wavefunction satisfying these relations exists, but does not allow us to calculate its form directly.One possible approach for modeling it is to use the hydrogenic wavefunctions and adjust their parameters such that they would satisfy the imposed overlap integrals (as employed in [25]).These modified hydrogenic wavefunctions, however, differ significantly from those of real carbon atoms bound within the graphene lattice which is limiting their predictive powers and makes the atomic contribution to the total electronic momenta unreliable.
Another approach is to use the Roothaan-Hartree-Fock wavefunctions fit to describe unbound carbon atoms in-stead of hydrogenic wavefunctions as the atomic basis.While this atomic basis does describe the individual carbon atoms better than the hydrogenic wavefunctions, they do not satisfy the overlap integrals given by the structure of the graphene lattice, creating an inconsistency in the implementation of the theory.In order to satisfy these relations, one would have to significantly distort the shape of the atomic orbitals, spoiling the original fit describing the carbon atom.
The problem underlying the form of the atomic orbitals within TB, together with the robust predictive powers of the ground-state electron density of DFT have led us to recommend DFT as the framework of choice for graphene-like DM detector modeling.We will further expand this topic and use DFT to obtain predictions for various possible detector setups and DM candidates in the associated Paper II.
The research software Darphene and an updated version of QEdark-EFT used to obtain the TB and DFT results respectively, will be made publicly available [39,40].
where v ⊥ el = v − q 2µχe − me with v being the DM initial velocity in the detector rest frame and µ χe the DM-electron reduced mass.c i 's are the effective couplings, and j χ is the DM spin which we typically set to 1/2. and Furthermore, to rewrite the above equations one can use the following relations where ∆E e is the energy transferred to the target electron.
Appendix B: The lattice structure of graphene Graphene is a two-dimensional hexagonal lattice of carbon atoms as illustrated in the left panel of Fig. 8.The distance between two neighboring carbon atoms is The two lattice vectors, which can also be seen in the left panel of fig.8 are such that The same figure also shows the vectors N i , pointing to a carbon atom's three closest neighbor atoms, The lattice vectors b i of the reciprocal lattice, illustrated in the right panel of Fig. 8, are The figure also shows the three high-symmetry point in the first Brillouin zone (BZ) (shaded in blue) in k-space, Appendix C: The tight-binding approximation

General review
In this section, we review the tight-binding approximation.In particular, we emphasize how the energy dispersion E i (k), appearing in Eq. ( 28), and wave function coefficients C ij (k) of Eq. ( 34) are evaluated in general using this technique [70].The translational symmetry of a lattice should also be reflected in the wave functions.In particular, the electron wave function Ψ(x) has to satisfy Bloch's theorem, Here, we introduced the translation operator T a along one of the lattice vectors a, and the lattice momentum k.
One way to write down a generic wave function satisfying Bloch's theorem can be obtained using the tight binding approximation.A tight-binding Bloch function Φ jk (x) is an approximation to the system's wave functions which is defined by summing up the wave functions of the j th atomic orbital of isolated atoms at their respective lattice site, This way, Bloch functions sum up the wave functions ϕ j (x) of N unit cells weighted with a lattice site dependent phase which ensures that Eq. (C1) is satisfied.Note that, even if the isolated atomic wave functions ϕ j (x) are normalized, the overlaps between neighboring wave functions generally render the Bloch functions as non-normalized.The actual electron wave functions of the material are linear combinations of the Bloch functions, mixing the different atomic orbitals (but not different lattice momenta), where the constant N k ensures that the wave function Ψ ik (x) is normalized.Using Schrödinger's equation, the energy values of these n states are Substituting Eq. (C3), the energy eigenvalues can be expressed as Here, we defined the coefficient vectors as well as the transfer integral matrix H(k), and the overlap integral matrix S(k).Their components are defined as As already seen in Eq. (C5), the normalization of Ψ ik (x) can be written in terms of the coefficient vector C i (k) and the overlap matrix S(k), The entries of the coefficient vector C i (k) are obtained using the variational principle by minimizing the energy eigenvalues E i (k), i.e.
This equation is equivalent to the following general eigenvalue problem, Non-vanishing eigenvectors can only be found, if the secular equation applies, For a fixed k, this polynomial of degree n can be solved for the n eigenvalues, i.e. the energy dispersion E i (k).Furthermore, the eigenvectors determine the electron wave functions Ψ ik (x).
Finally, the entries of the transfer integral matrix H(k), and the overlap integral matrix S(k) are often fixed to specific values, which ensure that the correct band structure of the material of interest is reproduced.These values are obtained either experimentally or from first-principles calculations.

Tight-binding approximation for graphene
We will apply the results of the previous section to the case of graphene.As illustrated in Fig. 8, the unit cell of graphene consists of two carbon atoms, denoted A and B. The relevant atomic orbitals of the carbon atoms on these locations are 2s, 2p x ,2p y , and 2p z (the 1s orbitals form a low energy core state and do not contribute to the valence band levels).In graphene, the first three orbitals combine or hybridize and form the so-called σ-bands and the 2p z orbitals hybridize to form the π-bands.As a result, there are n = 2(6) π(σ)-bonding atomic orbitals in each unit cell, and Eq.(C12) constitutes a two-(six-) dimensional eigenvalue problem for the π(σ) band.As we will describe in detail, the energy dispersion and wave function for the π electrons can be expressed analytically, whereas the σ electrons require numerical methods to solve the six-dimensional secular equation.

The π-electrons
The π-electrons are a hybridization of the atomic 2p z orbitals of carbon.In order to compute the corresponding energy bands by solving Eq. (C12), the main step is to evaluate the transfer integral and overlap integral matrices H(k) and S(k), which we defined in Eqs.(C7) and (C8).In this case, they are 2 × 2 matrices.Starting with the former, we substitute Eq. (C2) into Eq.(C7).The diagonal entries read If we only take the terms into account for which R k = R k and neglect sub-dominant contributions with R k = R k , we find By analogy, H BB = ε 2p .For the off-diagonal components, we can do a similar (nearest-neighbor) approximation.
Next, we only involve the contributions of the three nearest neighbors of each atom A, where the three vectors N k are given by Eq. (B3).We define the parameter t which is identical for all k due to the rotational symmetry of the 2p z wave function.This leaves us with The other off-diagonal is simply H BA = H * AB .In summary, the full transfer integral matrix for the π-electrons in the nearest-neighbor approximation is given by The function f (k) is defined as The square of this function can be evaluated as (C20) For the overlap matrix, the previous steps can essentially be repeated to find where Energy bands of π-electrons Given the explicit form of the two matrices, we can show that the corresponding eigenvalues solving the secular equation in Eq. (C12) are given by This energy dispersion is visualized in Fig. 9 for k ∈ BZ, as well as along the path between the high-symmetry points given by Eq. (B5), which is also indicated in the figure.In order to reproduce the band structure of graphene to a good degree, we used the parameters ε 2pz = 0 (by convention), s = 0.129 and t = −3.033eV[70].Since t is negative, the '+ solution is lower energy and corresponds to the bonding or valence π-band, whereas the '− solution is the anti-bonding or conduction π *band.The bands are degenerate at the high-symmetry point K.
Next, we turn our attention towards the C i (k) coefficients.The normalized eigenvectors corresponding to the eigenvalues of Eq. (C22) are obtained by solving Eq. (C11).For the π-electrons we find with Hence, the π-electron wave functions can be written as As opposed to the treatment in [25], we do not perform a nearest-neighbor approximation on this level, as it is not well-defined here.Instead it is necessary to sum over all N unit cells.However, the nearest-neighbor approximation can be applied when computing the norm of Ψ πk (x), Next, we perform the sum over k using the nearestneighbor approximation (but also taking next-nearest neighbors into account).In addition to the neighboring contributions, the first and last lines contain terms with R k = R k , each contributing with N cell to the final sum.Hence, we find Here, s denotes the overlap integral of the atomic orbitals at next-to-nearest neighboring sites.We find that the two leading terms are consistent with the general expression of Eq. (C9), and hence Next, we shift our attention from position space to momentum space.The Fourier-transformed Bloch wave functions are where is the conjugate momentum to x.Therefore, the π-electrons' wave functions in momentum space read (C32) For the evaluation of the exponential sum, we can follow the steps outlined in Sec.II D 1 and use the identity in Eq. ( 36).In the end, we find The σ-electrons The σ-electrons are in a superposition of the carbon atoms' 2s, 2p x , and 2p y orbitals.Hence, the transfer and overlap matrices are 6 × 6 matrices, which we can express in terms of four 3 × 3 sub-matrices, The diagonal sub-matrices are diagonal, with ε 2s = −8.87eVand ε 2p = 0.The off-diagonal matrices are given by with entries The other off-diagonal sub-matrix is given by the conjugate matrix, i.e.The energy bands of the σ-electrons are obtained by solving the secular equation i.e.Eq. (C12).Since we are dealing with 6 × 6 matrices, we use the numerical functionality of the Eigen library [50].The resulting energy bands for the σ-electrons are depicted in Fig. 10.The same Eigen function that computes the eigenvalues of the matrices in Eq. (C34), i.e. the energy bands of the σelectrons, by solving the secular equation also result in the six-dimensional eigenvectors C σi (k) and therefore the wave functions according to Eq. (C3).
The six Bloch wave functions are given by Finally, the Fourier transform of the wave function is given by Here, we again used R B j = R A j + δ.Just like in the case of the π-electrons, we use Eq. ( 36) to express the exponential sum in terms of a sum over the reciprocal lattice vectors G, Appendix D: Atomic wavefunctions The evaluation of the graphene response function requires a specific form of the atomic wavefunctions of carbon ϕ i (x) (and its Fourier transform ϕ i ( )).We present results for two particular choices.As proposed by Hochberg et al. [25], we start by approximating the wavefunctions of carbon with hydrogenic wave functions with a re-scaled Z eff factor.We improve upon this choice by using Roothaan-Hartree-Fock (RHF) wavefunctions for the ground states of carbon [51].In this appendix, we summarize the explicit wave functions both in position and momentum space and present a comparison.
The wave function of the atomic state (n, , m) in position space is given by where Y m l (x) are spherical harmonics, and R nl (r) is the radial component of the wave function.
Regardless of the choice of form for the radial component, the corresponding Fourier-transformed wave function in momentum space ϕ i ( ) for a given position space wave function ϕ i (x) is defined as which fixes the normalization of the wave function to For the evaluation of the graphene response function using the TB approximation, the required wavefunctions are those of the atomic orbitals 2s and 2p in the environment of a carbon atom.Additionally for the 2p orbitals, the crystal structure of graphene gives rise to the 2p x , 2p y , and 2p z orbitals.The corresponding wavefunctions are given by Y i (x) ≡ 3 4π Similarly, the relevant momentum wave functions can be written as and hence for the 2p i states we find ϕ 2pi ( ) = χ 2p ( )Y i ( ˆ ) .(D6)

Hydrogenic wave functions
We list the hydrogenic wavefunctions proposed in [25] to approximate the groundstate wave functions of carbon atoms.In position space, the wave functions are given by ϕ 2s (x) = (Z 2s eff ) 3 Using these expressions, we can evaluate all relevant atomic orbitals for the graphene response function by applying Eq. (D6).
When we compare the hydrogenic to the RHF wavefunctions in Fig. 11, we find the hydrogenic wavefunctions to be a poor match to the RHF wavefunctions, which are accepted as a good description for the atomic carbon groundstate wave functions.Instead, we conclude the necessity to use actual carbon atomic wavefunctions.This conclusion becomes even more robust when comparing the graphene response functions computed with TB and DFT.The first study of graphene targets for sub-GeV DM searches was published by Hochberg et al. [25].Therein, the authors chose a semi-analytic approach to describe the electron wavefunctions in graphene based on the tight-binding (TB) approximation.In reproducing their results, we noticed a number of deviations to our TB treatment of the electron wavefunctions in graphene.
• The modeling of the Bloch wavefunctions in the above-mentioned work does not satisfy the Bloch's theorem given in Eq. (C1).This also gives rise to a different normalization factor N k .
• While the hydrogenic wavefunctions proposed in [25] give overlap integrals in agreement with the TB theory, they differ significantly from more accurate RHF atomic wavefunctions of carbon atoms, as shown in App.D.
• After a careful evaluation of the formula for the DM induced electron ejection rate, we find an extra factor of 1/2 coming from normalizing to the number of unit cells instead of the number of carbon atoms in the system.
In this appendix, we review how we improved the TB treatment of graphene wavefunctions and how the updated electron ejection rate compares to the one presented by Hochberg et al. 4

Bloch states and their normalization
The "Bloch states" proposed in [25] are given by Compared to the expression of Eq. (C2), these states, while capturing the nearest-neighbour approximation, do not formally correspond to Bloch wavefunctions.Consequently, they do not satisfy Bloch's theorem or rigorously describe a periodic system, see Eq. (C1).One consequence of this choice of Bloch states is a deviating normalization factor N k .As an example, using the Bloch states by Hochberg et al., it is possible to compute the exact normalization factor for the π-electrons, While this factor correctly normalizes the electron wavefunctions involving the Bloch states of (E1), it deviates from our respective expression for the π-electrons in Eq. (C29).Furthermore, our normalization constant is consistent with Eq. (C9), i.e. the general expression for the normalization constant for Bloch wavefunctions as given by Eq. (C2).Similar arguments hold for the σ-electrons.[25]).The solid lines show the total spectrum, whereas the dashed (dotted) lines show the contributions of the π-(σ-) electrons.For this figure, we used the SHM parameters from [25] in order to facilitate the comparison (Hochberg et al. use v0 = 220 km/s as opposed to our choice of v0 = 238 km/s.)

Atomic wavefunctions
In [25], Hochberg et al. propose to describe the atomic wavefunctions of carbon by using hydrogenic wavefunctions with a re-scaled Z eff factor.The re-scaling ensures that the overlap integrals of the wavefunctions are consistent with the TB parameters that reproduce the band structure of graphene listed in Tab.II.For completeness, we list these wavefunctions in App.D 1.
In contrast, we model the carbon wavefunctions using Roothaan-Hartree-Fock (RHF) wavefunctions [51], which we summarize in App.D 2. By comparison, we find that the re-scaled hydrogenic wavefunctions are a poor approximation for the required ground state wavefunctions of atomic electrons in carbon, as can be seen in Fig. 11.
While the RHF wavefunctions are a better description of carbon electrons, it should be noted that their overlap integrals are not consistent with the overlap parameters for graphene in Tab.II.A re-scaling of the RHF wavefunctions similar to the approach by Hochberg et al. is possible, but spoils the accurate description of electrons in atomic carbon.However, this seems to be a general feature of evaluating electron wavefunctions in the TB approximation.

Electron ejection rate
In [25], the total rate of DM-induced electron ejections in graphene is given as where Here, we use the notation of Hochberg et al.This needs to be compared to our Eqs.( 31) and (32).Here, we use the replacement Eq. ( 44) for two-dimensional targets, and further identify to facilitate the comparison.We find agreement between our expressions with one exception.Instead of the number of carbon atoms N C , we find that the electron ejection rate is proportional to the number of unit cells N cell .Hence our expressions for the electron ejection rates differ by a factor of 2. In Fig. 12, we compare the energy spectrum for a DM particle of 100 MeV mass using the TB approach presented by Hochberg et al. and compare it to the improved version presented in this work.Note that the interaction model used by Hochberg et al. corresponds to O 1 interactions in our general framework.We find that our predicted spectrum is about one order of magnitude lower than the one presented in [25].

FIG. 1 .
FIG.1.The partially integrated graphene response function evaluated with TB (left) and DFT (right) as a function of x and y .In the left panel, we also show the contributions of each of the electron bands.We set z = 91 eV, such that the vector lies almost in the plane of the graphene sheet.The stripe-like structure of the DFT response is an artifact of the grid sampling in the reciprocal space and is integrated out when observables are evaluated.

FIG. 2 .
FIG. 2. W ( , Ee) integrated between Emin < 0 and 0 and over the azimuthal angle of , φ for TB (left) and DFT (right).The first row is for Emin = −5 eV, including only electrons accessible to low mass DM.The second row has Emin = −20 eV and also includes electrons accessible to heavier DM.As in Fig.1, the distorting effect of the finite sampling of the reciprocal space can be seen for the case of DFT.When evaluating the final state observables, these distortions are washed out.

FIG. 3 .
FIG.3.Graphene response function and its dependence on the initial state momentum integrated over various regions of initial state energy and all directions.This plot shows how much different electron energies (that are accessible to different DM candidate masses) contribute to the momentum distribution of the target for TB (left) and DFT (right) simulations.

FIG. 4 .
FIG.4.Comparison of the partially integrated response function as a function of the initial electron momentum, , between DFT (blue) and TB(red).The solid lines illustrate the response function for momenta perpendicular to the graphene sheet.The dashed lines show the dependence of W on momenta parallel to the graphene sheet, where we average over the in-plane directions.The lower panel shows the ratio.As indicated by the gray band, for momenta below ∼ 15 keV, the two predictions are within a factor of 2. Above, the TB approximation predicts a significantly higher response function.

FIG. 8 .
FIG.8.The honeycomb lattice of graphene (left) and the reciprocal lattice (right), together with their respective lattice vectors ai and bi.The red (blue) shaded area is the unit cell (Brillouin zone) of the (reciprocal) lattice.

FIG. 10 .
FIG.10.Energy bands of the π and σ-electrons obtained using the tight-binding approximation.Solid lines show the valence bands, dashed lines are the conduction bands.
Appendix E: Comparison of our TB treatment toHochberg et al. (2017)

m 1 ]FIG. 12 .
FIG.12.Comparison of the energy spectrum between our TB results and the results by Hochberg et al. (spectrum digitized from Fig.1of[25]).The solid lines show the total spectrum, whereas the dashed (dotted) lines show the contributions of the π-(σ-) electrons.For this figure, we used the SHM parameters from[25] in order to facilitate the comparison (Hochberg et al. use v0 = 220 km/s as opposed to our choice of v0 = 238 km/s.) . 7. Daily modulation pattern for graphene sheets obtained with DFT.The colors in each panel correspond to the interaction types indicated in the legend, and the top left, top center, top right, bottom left, bottom center and bottom right corresponding to DM masses of 2 MeV, 5 MeV, 10 MeV, 20 MeV, 50 MeV and 100 MeV, respectively.Note that the y-axis differs between the top and bottom plots.Note that for DM masses around 10 MeV, the modulation pattern is similar for all the considered interactions, indicating that the expected result is largely model-independent. FIG

TABLE II .
S BA = S † AB .The off-diagonal sub-matrices H AB and H BA are obtained by replacing S → H and S → H.The numerical parameters are given in Table II.Parameters of the overlap and transfer matrices for the σ-electrons