Atomistic Band-Structure Computation for Investigating Coulomb Dephasing and Impurity Scattering Rates of Electrons in Graphene

In this paper, by introducing a generalized quantum-kinetic model which is coupled self-consistently with Maxwell and Boltzmann transport equations, we elucidate the significance of using input from first-principles band-structure computations for an accurate description of ultra-fast dephasing and scattering dynamics of electrons in graphene. In particular, we start with the tight-binding model (TBM) for calculating band structures of solid covalent crystals based on localized Wannier orbital functions, where the employed hopping integrals in TBM have been parameterized for various covalent bonds. After that, the general TBM formalism has been applied to graphene to obtain both band structures and wave functions of electrons beyond the regime of effective low-energy theory. As a specific example, these calculated eigenvalues and eigen vectors have been further utilized to compute the Bloch-function form factors and intrinsic Coulomb diagonal-dephasing rates for induced optical coherence of electron-hole pairs in spectral and polarization functions, as well as the energy-relaxation time from extrinsic impurity scattering of electrons for non-equilibrium occupation in band transport.


Introduction
Very recently, a generalized parameter-free quantum-kinetic model [1,2] based on many-body theory [3,4] has been developed, which is self-consistently coupled with Maxwell equations [5] for an interacting electromagnetic field and with Boltzmann transport equation [6] for a conduction current, as illustrated in Figure 1. Here, being an off-diagonal element in a density matrix, the induced quantum coherence for electron-hole pairs leads to a macroscopic optical polarization field [1] included in the Maxwell equations. Meanwhile, the modified electric field determined from the Maxwell equations can also change the microscopic quantum coherence [1] of electron-hole pairs. In this way, a selfconsistent loop is constructed between electrons in the quantum-kinetic model and electric field in the Maxwell equations. This theory aims at enabling first-principles computations of ultra-fast dynamics for non-thermal photo-generated electron-hole pairs in undoped semiconductors [1,2]. At the same time, this a theory is also able to simultaneously describe electromagnetic, optical and electrical properties of crystal materials and their interplay all together. More importantly, the numerical output of this first-principles dynamics model can be utilized as an input for material optical and transport properties to be fed into a next-stage simulation software facilitated by finite-element methods, such as COMSOL The first-principles computation of electron Bloch wave function and band dispersion of a targeted material can be performed by employing the well-known Kohn-Sham densityfunctional theory [10]. Meanwhile, the tight-binding model [11][12][13][14][15] for solid crystals is usually considered as an alternative approach for computing electronic band structure using an approximate set of orbital wave functions based upon superposition of bondorbital states for isolated atoms sitting at different lattice sites. In fact, this method is closely related to the linear combination of atomic orbitals method [16] adopted commonly in quantum chemistry. Such a real-space tight-binding model can be applied to a lot of solids, even including a magnetic field, Ref. [17] and it is proved giving rise to good qualitative results [18]. Moreover, this method can be combined with other models to produce better results whenever the tight-binding model fails. Here, we would like to emphasize that although the tight-binding model is only a one-electron model in nature, it indeed provides a basis for more advanced computations [11], such as the computation of surface states, application to various kinds of many-body problems, and quasi-particle calculation [19].
Historically, the family of carbon-based materials can be characterized into two distinct crystal forms, i.e., the isotropic diamond and anisotropic graphite. Recently, their allotropes, such as fullerenes and carbon nanotubes, entered into play and expanded to graphene, which is a unique material consisting of a two-dimensional lattice of carbon atoms with a honeycomb symmetry. Graphene stands for an physically interesting system [20,21], and becomes very promising for future device applications. On the atomic level, e.g., densityfunctional theory, electron certainly follows the Schrodinger equation. However, by using an approximate effective-mass Hamiltonian [22,23] for low-energy electrons near the K or K valley, the quasi-particles are found to satisfy the relativistic Dirac equation for massless fermions. Today, the extensive investigations on various graphene systems have turned into a broad research field for qualitatively new two-dimensional systems [24]. Up to now, the basic properties of novel 2D allotropes of carbon, including graphene [22,23], graphene bilayer [25][26][27], multi-layer graphene [28,29], graphene on a silicon carbide substrate [30], are well known and the basis of graphene physics becomes well established.
The rest of paper is organized as follows. In Section 2, we present a general description of tight-binding model for novel two-dimensional materials. Section 3 is devoted to discuss the Slater-Koster approximation for bonding parameters and bonding integrals. We acquire the parameter values in Section 4 and obtain graphene wave functions and band structures. We study the Coulomb diagonal-dephasing rate of electron-hole pairs in undoped graphene in Section 5, as well as the impurity scattering rate of conduction electrons in Section 6, respectively. Finally, a brief summary is presented in Section 7 along with some remarks.

General Description of Tight-Binding Model
For completeness, we start with tight-binding model [14] for computing complete band structures of two-dimensional materials. The advantage of tight-binding model is easily incorporating a magnetic field through the so-called Peierls substitution in the phase of a hopping integral [51]. In quantum mechanics, the single-electron static Schrödinger equation is written as [52] where Φ k (r) is the Bloch wave function, E k the eigen-energy, and k is the wave vector of electrons within the first Brillouin zone of two-dimensional materials. The Hamiltonian operatorĤ in Equation (1) in which the kinetic-energy operatorĤ 0 iŝ with free-electron mass m e , while the potential energy V(r) for an electron within the lattice of two-dimensional materials is given by [11] V(r) = V ion (r) + ∆V L (r) (4) with V ion (r) and ∆V L (r) specifying the potentials of a single ion and that for the rest of ions, respectively. The Bloch wave function Φ k (r) of electrons in Equation (1) can be decomposed into a linear combination of a set of orbital wave functions φ β,k (r) within the first Brillouin zone, leading to [11] where the index β labels all the atomic orbitals of the lattice of two-dimensional materials.
The expansion coefficient C β introduced in Equation (5) can be decided from where the orthonormal property for the set of orbital wave functions φ β,k (r) has been adopted. Applying the method of linear combination of atomic orbitals (LCAO) of all ions on the lattice [16], we further express each orbital wave function φ β,k (r) in Equation (5) by a linear combination of bond-orbital states ψ β (r − R j ) within a unit cell in real space, namely where j is the index for all bonded lattice ions, R j the lattice-ion position vector, and N the total number of atoms within the unit cell. Here, is termed as the localized Wannier function for the β orbital of a bonded lattice ion at the site R j , which satisfies the single-ion Schrödinger equation [16] with ε j,β being the βth energy levels of electrons within an ion at the lattice site R j . Combining results in Equations (5) and (7), we acquire the following full LCAO expansion of a Bloch wave function [11] with C β; j (k) = C β exp(ik · R j ). At the same time, using Equation (5), we find from or equivalently, the following eigenvalue equation As a result, the eigenvalue E n (k) can be determined from the secular determinant of Equation (12) for any given k, yielding and the orthonormal-eigenvectors C n,β are also obtained, corresponding to the eigenvalue E n (k) at given k, where the index n labels different quantized energy bands of two-dimensional materials. Explicitly, using Equation (7), we obtain the Hamiltonian matrix elements in Equation (12) as [11] H αβ (k) in which R ij = R j − R i , and In fact, we know from Equation (9) that where (ε j,β + C Σ ) represents the site energy, and t αβ (R ij ) is usually called the two-center (or hopping) integral [14]. As a final step, with the help from Equation (10), we arrive at the full expression for Hamiltonian matrix elements, given by where the primed summation in the second term of the right-hand side of the last equation excludes the contribution from j = j , and C n,β can be obtained from the calculated eigenvector from Equation (12). The matrix elements for other physical operators can be computed in a similar way.

Slater-Koster Approximation for Hopping Integrals
To seek for the feasibility of fast numerical computation, we introduce a parameterized process for the tight-binding model described in Section 2. For the Coulomb interaction between electron and ion within an atom, the potential field presents a spherical symmetry. Therefore, the energy levels labeled by the radial quantum number n = 1, 2, · · · will degenerate with the angular-momentum quantum number = 0, 1, · · · , n − 1, as well as the magnetic quantum number m = − , · · · , 0, · · · , [52]. Consequently, there exists a total orbital degeneracy n 2 (excluding the spin-degeneracy). Customarily, we specify these orbitals by = 0, 1, 2, 3, · · · for {s, p, d, f , · · · } orbitals.
In order to describe the chemical bonds between a pair of atoms inside a lattice, we often adopt the concept of overlapping electronic orbitals {s, p, d, f , · · · }. To further specify the spatial direction of the chemical bonding between two atoms at the lattice sites R i and R j , we have to rely on three directional cosines , m, n, as defined in Figure 2.
Considering s and p orbitals as an example, we display their possible bonding potentials V , ;σ(π) in Figure 3 for s, p orbitals and four different configurations, including σ and π bonds. Meanwhile, we also list six different π, σ, δ bonding configurations in Figure 4 for s, p, d orbitals.
To speed up numerical computations, the bonding potentials V , ;σ(π) for , = s, p, d in Figures 3 and 4 are usually parameterized as: [53] where d and r d represent the bonding length and atomic radius, and ξ = σ, π, δ are for various bond configurations. Here, the dimensionless bonding parameters η , ;ξ for different bonding types are listed in Table 1.
By using these parameterized bonding potentials V , ;ξ , V ,d;ξ and V d,d;ξ , we are able to compute further the hopping integrals t αβ (R ij ) based on the Slater-Koster approximation [14], and some commonly-used results are shown in Table 2.

Tight-Binding Model for Graphene Band Structure
To seek for an application, we use the general theory, as developed in Sections 2 and 3, for novel two-dimensional graphene material in order to obtain its electronic wave functions and band structures for the full first Brillouin zone [54]. In this way, we are able to study scattering dynamics with respect to high-energy electrons in graphene resulted from Coulomb interactions between either pair of electrons or between electrons and ionized impurity atoms.
Monolayer graphene displays a hexagonal (or honeycomb) lattice structure of carbon atoms, as illustrated in Figure 5, where each carbon atom is connected by σ covalent bonds with its three nearest neighbors. The electronic orbitals of a carbon atom are characterized as 1s 2 2s 2 2p 2 . However, the unique energy difference between the 2s and 2p orbitals favors the appearance of a mixed state of these two orbitals. The first-principles density-functional calculations reveal that it becomes energetically favorable to move an electron from the 2s orbital to the 2p orbital in this mixed state. Since the 2p orbitals include 2p x , 2p y , 2p z , as a result, each of these three 2p orbitals will accommodate one electron, leading to the x-y orbitals within the plane of the lattice, as well as the z orbital out of the lattice plane. Here, two electrons in the mixed x-y orbitals form the higher-energy σ bonds, while the remaining electron in the z orbital leads to the lower-energy π bonds, i.e., a side-on overlap of the 2p-orbital wave functions. Consequently, these π-bond electrons give rise to the low-energy bands of graphene and will be studied exclusively based on a tight-binding model. From Equation (7), we know the wave function for π-bond (p z -orbital) electrons in graphene can be expressed as where 2 represents the Bravis lattice-site vectors as indicated in Figure 5, and indexes A, B refer to two sublattices of graphene. By including both sublattices A and B, we have where a k and b k are two elements of the eigenvector corresponding to the eigenvalue equation with respect to two sublattices. Specifically, from Equations (11) and (20), we arrive at the matrix-form Schrödinger equation represents the eigen-energies of π-bond electrons with n = 1, 2 labeling two graphene lowenergy bands determined by the secular determinant: As in Equation (7), we can rewrite the orbital wave function φ p z ,k (r) in Equation (20) approximately only by its near-neighbor decomposition, yielding and then, the eigenvalue equation where N c represents the number of near-neighbor atoms within a unit cell, ∆ j stands for the lattice vectors of the near-neighbor atoms relative to the sublattice site R , and a jk = a k exp(−ik · ∆ j ). Moreover, we find H jj (k) = ε 2 S jj (k) + t jj (k), where ε 2 stands for the second energy level of electrons within a carbon atom, is the overlap integral, while is the hopping integral. For simplicity, we would omit the orbital index p z from now on. Without loss of generality, we can assume that the vectors that connect sublattice A site to the equivalent site on the B sublattice is δ 3 , as seen in Figure 5. As a result, the hopping and overlap amplitudes between the nearest neighbor (nn) and the next-nearest neighbor (nnn) can be computed explicitly from Equations (23) and (24), leading to where a 3 ≡ a 1 − a 2 , γ(k) = 1 + exp(ik · a 1 ) + exp(ik · a 2 ), and the hopping and overlap integrals are calculated as Particularly, the results for these tight-binding model parameters in Equation (26) for band structures are presented in Table 3, which have been computed from listed bonding parameters in Table 1 and bonding integrals in Table 2.
Finally, from the eigenvalue equation Det{ t (k) − E(k) S (k) } 2×2 = 0 in Equation (21) for , = A, B, we obtain an explicit expression where, by setting s nnn = 0, we have three coefficients This leads to the explicit solution of Equation (27), namely [55] where λ = ±1 correspond to valence (−1) and conduction (+1) bands, respectively, and By using the result in Equations (29) and (30) can be rewritten as By setting C p + ε 2 = 0 as the reference point for energy, the result in Equation (31) is plotted in Figure 6 by employing the graphene structural parameters listed in Table 3. Furthermore, by using the result in Equation (31), two elements of the eigenvector, a λ k and b λ k , are found to be As known experimentally, both the nearest-neighbor (nn) overlap and the next-nearestneighbor (nnn) hopping integrals are much smaller than the nearest-neighbor (nn) hopping integral. By neglecting some constants, the dispersion in Equation (31) can be further simplified as where t nnn = t nnn − s nn t nn is the corrected hopping amplitude.

Coulomb Diagonal-Dephasing Rate for Optical Coherence in Undoped Graphene
The quantum coherence of electrons is associated with the off-diagonal elements of their density matrix. The presence of an external field can induce coherence between two quantum states of electrons if the field frequency matches the energy separation between the two relevant electronic states. Dephasing refers to a physics mechanism which recovers classical behavior from a quantum system, and it quantifies the time required for electrons to lose their field-induced quantum coherence. Diagonal-dephasing rate connects to the ways in which coherence caused by perturbation decays over time, and then, the system goes back to the state before perturbation [1]. This is an important effect in molecular and atomic spectroscopy, and also in condense-matter physics of mesoscopic devices.
In order to demonstrate the significance of band-structure computation with a tightbinding model on dynamical properties of electrons in graphene, we first study Coulomb diagonal-dephasing (CDD) rate for induced optical polarization of thermally-excited electrons and holes around the Dirac point in an intrinsic (or undoped) graphene sample. For undoped graphene, conduction electrons can be introduced by a photo-excitation process [8], giving rise to equal number of electrons and holes n e = n h ≡ n 0 , where n 0 represents the areal density of photo-excited carriers. For non-equilibrium photo-carriers under a transverse optical field, its induced optical coherence in steady states decays [1] with the sum of CDD rates ∆ e (k) and ∆ h (k) for electrons (e) and holes (h), respectively. These two rates determine the inhomogeneous line-shape of a resonant interband-absorption peak at hω = ε e k + ε h k for vertical transitions of electrons with their kinetic energies ε e,h k in valence and conduction bands.
As illustrated by Feynman diagrams [3] in Figure 7, the CDD rate ∆ e (k) of electrons is calculated as [1,44] ∆ e (k) = 8π where both spin and valley degeneracies are included, A represents the surface area of graphene sample, the first and second terms correspond to the left and right panels of Figure 7, and both scattering-in and scattering-out contributions [44] are taken into consideration in these two terms. Moreover, ε e,h k in Equation (34) where both spin and valley degeneracies are included and µ e = µ h in our case. Furthermore, in Equation (34), L 0 (a, b) = (b/π)/(a 2 + b 2 ) is the Lorentzian line-shape function, Γ e,h are inverse lifetime of unperturbed electrons or holes, and Γ eh = (Γ e + Γ h )/2.

Diagonal Dephasing Rate
coulomb coupling between pair of electrons in a scattering event for inverse electron lifetime coulomb coupling between an electron and a hole in a scattering event for inverse electron lifetime In addition, we have introduced in Equation (34), as well as in Equation (40) below, the Coulomb-interaction matrix elements, given by [56] V ee k,k 1 ; where u c (q) = e 2 /[2 0 r (q + q 0 )]) in Equation (36) is the two-dimensional Fourier transformed Coulomb potential ∼ 1/r including static screening, 0 represents the vacuum permittivity, and r = 2.4 is the average dielectric constant of the host material. Additionally, q 0 stands for the inverse Thomas-Fermi screening length, and can be given by a semi-classical model as [57] where both spin and valley degeneracies have been included.
Furthermore, the introduced F (s) k, k (q) in Equation (36) with s = c, v represents the Bloch-function form factor, calculated as [57] F (s) where the Bloch functions Φ c,v k (r) in Equations (5) and (22) have been employed. In Equation (38), N c represents the number of near-neighbor atoms within a unit cell, ∆ j stands for the lattice vectors of the near-neighbor atoms relative to the sublattice site R , and a (s) k are two column eigenvectors in Equation (32) for s = c, v. The Wannier-function structure factor W s (q, R + ∆ jj ) in Equation (38) is defined as where R = R − R and ∆ jj = ∆ j − ∆ j . In fact, Equations (34) and (36)- (39) are the key results in this paper for connecting the calculated tight-binding wave functions and band structures to a quantum-statistical theory for graphene optical properties. Similarly, as illustrated by Feynman diagrams [3] in Figure 8, the CDD rate ∆ h (k) of holes takes the form [1,44] ∆ h (k) = 8π coulomb coupling between pair of holes in a scattering event for inverse hole lifetime coulomb coupling between a hole and an electron in a scattering event for inverse hole lifetime Computationally, the π-electron band structure of graphite can be obtained by employing the nearest-neighbor tight-binding model [58,59]. For graphene, the reciprocal lattice in the wave-vector space also acquires the hexagonal symmetry, same as that in real lattice. Moreover, the low energy bands are found linear and isotropic near the corners of the first Brillouin zone or K point. Such K-point linear bands become essential for the low-energy (or small wave-number) excitation of electrons. The calculated energy dispersions by diagonalizing the 2 × 2 Hamiltonian matrix are given by [58,59] ε e,h k = ± where γ 0 = 2.4 eV is the hopping integral between the nearest-neighbor atoms, b = 1.42 Å is the C-C bond length, and signs ± represents conduction (+) and hole (−) bands, respectively. Meanwhile, the corresponding spinor-type Bloch wave functions are found to be where as shown in Equation (20), U k (r) and U (2) k (r) are two sublattice Bloch functions built from the superposition of the periodic 2p z orbitals, Ref. [59] and θ k = tan −1 (k y /k x ) is the angle between the wave vector k and x-axis. As in Equation (22), we can further express the 2p z atomic orbital by means of a generalized hydrogen-like wave function, given by [60] ψ p z (r) = C 0 r cos θ e −Z * r/2a 0 , where C 0 is a normalization factor, a 0 the Bohr radius, and an effective nucleus charge number Z * is 3.18.
In particular, the structure factor introduced in Equation (38) can be calculated explicitly as where s = c, v for Bloch wave function. Moreover, the Bloch-function structure factor in Equation (44) takes the form [60] φ (s) where tight-binding function ψ p z (r) is given by Equation (43), and the signs (±) correspond to conduction (+) and valence (−) bands, respectively [59]. For intrinsic graphene, we have chemical potential µ e = µ h = 0 [61]. However, there is still a finite intrinsic areal density n i ≈ (π/6) (k B T/hv F ) 2 due to thermal excitation of electrons and holes at finite temperatures T. In fact, we find f e k = f h k = 1/2 at the K valley or k = 0. Here, the calculated CDD rates from Equations (34) and (40), respectively, for electrons ∆ e (k) and holes ∆ h (k) are presented in Figure 9a at T = 77 K and in Figure 9b at T = 300 K. Since f e,h k ∼ exp(−ε e,h k /k B T) as ε e,h k k B T, the thermal occupations of electron and hole states will be limited mostly to wave numbers close to the K valley due to their lower kinetic energies ε e,h k around k = 0, as seen in Figure 6. The Coulomb diagonal-dephasing rates ∆ e,h (k) presented in Figure 9a,b quantifies an amplitude-decay process of induced electron-hole optical coherence with wave vector k by an optical field towards the state before external perturbation. Furthermore, the Coulomb off-diagonal-dephasing rates Λ e,h (k, q) reveals deformations of induced opticalpolarization waves with different wave vectors k + q [8].
Considering the fact that major occupations of electrons and holes are accumulated around k = 0, we have f e,h k ≈ 0 only if k is large. As a result, we find from Equation (34) that f e k 1 −q f e k+q (1 − f e k 1 ) 1 at k = 0 since we require 1 − f e k 1 ≈ 1 for large k 1 , f e k+q ≈ 1 for small q, and f e k 1 −q ≈ 1 for both large q and k 1 , which, however, cannot be satisfied simultaneously. Similar conclusion can also be drawn for the second term in Equation (34), where we find f e Combining these two facts together, we expect that a dip will occur at k = 0 for the Coulomb diagonal-dephasing rate ∆ e (k), as seen in Figure 9a. Moreover, the observed anisotropic energy dispersion in Figure 6a along the K-M and K-Γ directions directly leads to a staircase-like feature in Figure 9a for both ∆ e (k) and ∆ h (k). As temperature T is raised from 77 K in Figure 9a to 300 K in Figure 9b, the thermally-excited areal densities of electrons and holes are increased with T 2 ; therefore, the Coulomb interaction (∝ T 4 ) between electrons and holes, as well as the Coulomb interaction among electrons or holes, will be enhanced greatly. Consequently, we find that both ∆ e (k) and ∆ h (k) are enhanced by a factor of 2.4, in addition to amplified depth of the dip at k = 0. Furthermore, different structural factors in Equation (39), corresponding to ± signs for conduction and valence bands, give rise to a slightly larger value of ∆ h (k) in comparison with that of ∆ e (k), as well as different dispersion features around the K valley for ∆ e (k) and ∆ h (k). These two computed Coulomb diagonal-dephasing rates can be physically applied to the spectral [32] and polarization [1,36] functions in order to study transport and optical properties of graphene material.

Carrier Energy-Relaxation Rate in Doped Graphene
In condensed-matter physics, the microscopic energy-relaxation time usually refers to a measure of the time it requires for one electron in the system to be significantly affected by the presence of other electrons, lattice vibrations, and randomly-distributed ionized impurity atoms in the system through an either scattering-in or scattering-out process mediated by electron-electron, electron-phonon and electron-impurity interactions, respectively. Since the microscopic energy-relaxation time is assigned to a specific electronic state, we are able to define a thermally-averaged energy-relaxation time through the diagonal density-matrix elements of electrons for all electronic states. In this way, one can reveal unique temperature dependence of this macroscopic energy-relaxation time and utilize it for simplifying the well-known Boltzmann transport equation within the relaxation-time approximation [44].
By going beyond the intrinsic graphene samples, we would like to investigate further the impurity scattering of electrons in extrinsic (or doped) graphene materials. In parallel with the discussion on scattering rates in Section 5, we present here the calculations for intraband-scattering of electrons by randomly-distributed impurities. Results for intrabandscattering of holes can be obtained in a similar way.
By using the detailed-balance condition, the microscopic energy-relaxation time τ rel (k) of electrons in the presence of randomly-distributed ionized impurities can be calculated according to [44] 1 where the scattering-in rate for electrons in the final |k state is whereas the scattering-out rate for electrons in the initial |k state takes the form Here, n im represents the areal density of ionized impurity atoms in the crystal, and |U im k, k (q)| 2 comes from the randomly-impurity scattering of electron in the second-order Born approximation [48,62]. Explicitly, the random impurity-interaction matrix elements are calculated as where Z * is the charge number of ionized impurity atoms. Substituting Equation (49) back into Equation (47), we obtain Using the inverse microscopic energy-relaxation time in Equation (46), we can further calculate the macroscopic thermally-averaged energy-relaxation time τ rel (T) as a function of temperature T, yielding [44] 1 τ rel (T) Actually, the results in Equation (46) and in Equations (50)-(52) demonstrate the approach for relating the computed tight-binding wave functions and band structures to graphene transport properties described by a many-body scattering theory. This calculated relaxation time in Equation (52) can be employed for building up different orders of moment equations [63] based on semi-classical Boltzmann transport equation [6] under the relaxation-time approximation [44]. Here, the zeroth-order moment equation [63] grantees the conservation of conduction electrons and allows us to find the chemical potential of electrons, as in Equation (35), for given areal doping density and temperature. Moreover, the first-order moment equation [63] makes it possible to find transport mobility and conductivity [64] for bias-field driven conduction electrons.
For doped graphene, we have Fermi energy E F =hv F √ πn 0 at low temperatures, Ref. [61] where n 0 represents the areal electron density from doping, i.e., n 0 = n im for completely ionized doping atoms. For low temperatures with k B T E F , we have is a unity step function and k F = √ πn 0 is the Fermi wave number.
Physically, the Coulomb diagonal-dephasing rates Γ(k) = ∆ e (k) + ∆ h (k) in Figure 9 describes a decay process of induced electron-hole optical coherence, which is induced by an optical field over time, towards the state before perturbation. On the other hand, the electron energy-relaxation rate 1/τ rel (T), determined by Equations (46) and (52), reflects the time, which is a quantum-statistical average over all occupied states of electrons, needed for recovering from a non-equilibrium-state occupation after an external perturbation to an initial thermal-equilibrium-state occupation before external perturbation via an elastic electron-impurity scattering process. Therefore, these two rates, as shown by Figures 9 and 10, respectively, represent two fundamentally different microscopic physics mechanisms. Figure 10. Calculated average energy-relaxation rate 1/τ rel (T) from Equation (52) as a function of temperature T due to elastic scattering of doped electrons with impurities in graphene material, where r = 2.4, Γ e = Γ h = 0.01 meV, Z * = 1, doped electron areal density n 0 = 1 × 10 11 cm −2 , and impurity areal density n im = n 0 are assumed.
As seen from Figure 10, we find the electron energy-relaxation rate 1/τ rel (T) reduces with increasing temperature T due to enhanced screening effect on Coulomb interaction u c (q) between two electrons or the rising of q 0 in Equation (37) with T, which implies that we have to wait a longer time τ rel (T) for our system returning to its initial thermalequilibrium state at an elevated temperature. Furthermore, using the second-order Boltzmann moment equation [44], we would emphasize that this average energy-relaxation time τ rel (T), as determined from Equations (46) and (52), is directly associated with the mobility of transport electrons limited by elastic scattering from existence of impurities in the system.

Conclusions and Remarks
In conclusion, by introducing a generalized first-principles quantum-kinetic model coupled self-consistently with Maxwell and Boltzmann transport equations, we demonstrate the importance to incorporate inputs from first-principles band-structure computations for accurately describing non-equilibrium optical and transport properties of electrons in graphene. Generally speaking, the physical properties of an active material in a device are determined by both underlined band structures of involved materials and non-equilibrium responses to various external impulses.
In this study, we initialize with the tight-binding model for investigating band structures of solid covalent crystals by means of localized Wannier orbital functions, and further parameterize the hopping integrals in the tight-binding model for different covalent bonds. After that, we apply the general tight-binding-model formalism to graphene in order to acquire both band structures and wave functions of electrons within the whole first Brillouin zone of two-dimensional materials. For illustrating their significance, we utilize them to explore the intrinsic electron-hole Coulomb diagonal-dephasing rates used for spectral and polarization functions of graphene materials, and meanwhile, the energy-relaxation rate from extrinsic elastic scattering by impurities for transport mobility of doped electrons in graphene.
Theoretically, our current theory is capable of first-principles calculations of ultra-fast dynamics for non-thermal photo-generated electron-hole pairs. Simultaneously, this a theory also enables to describe electromagnetic, optical and electrical properties of semiconductor materials all together, as well as their interplay. Technologically, in combination with first-principles band-structure computations, the numerical output of current firstprinciples dynamics model can be used as an input for material optical and transport properties and put into a next-step simulation software, such as COMSOL Multiphysics, for a target device. Consequently, device characteristics can be predicted accurately for numerical bottom-up design and engineering.  Data Availability Statement: The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.