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

Abstract

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.

1. 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 self-consistent 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 Multiphysics [7], for devices with various configurations. Consequently, device characteristics can be accurately predicted beyond the linear-response regime [3,4] for numerical bottom-up design and engineering. However, such a quantum-kinetic model itself requires an input from wave functions and band structures associated with different host materials in devices.

In Figure 1, we introduce the product of field frequency ($\omega$) with the carrier momentum-relaxation time (${\tau }_p$). The situations with $\omega {\tau }_p\gg 1$ and $\omega {\tau }_p\ll 1$ correspond separately to optical and bias field regimes, while $\omega {\tau }_p\approx 1$ uniquely specifies the terahertz regime with dual optical and bias field characteristics. The bridging connection between the Maxwell [5] and semiconductor Bloch [8,9] equations is provided by the induced optical-polarization field $\mathit{P}(\mathit{r},t)$ as a quantum-statistical average of the electric-dipole moment with the induced microscopic optical coherence ${\mathit{p}}_j(\mathit{k},t)$ with j the band index. The bridging connection between the Maxwell [5] and Boltzmann transport [6] equations, on the other hand, is fulfilled by the optically-induced magnetization field $\mathit{M}(\mathit{r},t)$ as a quantum-statistical average of the induced microscopic magnetic-dipole moment ${\mathit{m}}_j(\mathit{k},t)$ from spins or orbital angular momentum. Finally, the bridging connection between the semiconductor Bloch [8,9] and Boltzmann transport [6] equations is facilitated by the bias-induced macroscopic center-of-mass drift velocity ${\mathit{v}}_d\left(t\right)$ as a non-equilibrium quantum-statistical average of the microscopic electron group velocities ${\mathit{v}}_j\left(\mathit{k}\right)$ from multi-band dispersions for modifying optical-transition properties of driven carriers within the center-of-mass frame due to relative scattering motions of carriers.

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 density-functional 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 bond-orbital 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., density-functional 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\prime$ 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.

In recent years, by using the low-energy Dirac Hamiltonian [4], we have extensively explored varieties of dynamical properties of electrons in graphene and other two-dimensional materials, including Landau quantization [18,31,32,33,34,35], many-body optical effects [36,37,38,39,40,41], band and tunneling transports [42,43,44,45,46,47,48,49,50], etc. In this paper, we particularly focus on the application of computed electronic states and band structures from a tight-binding model to the calculations of Coulomb and impurity scatterings of electrons in graphene on the basis of a many-body theory [3,4], where the former and latter determine the lineshape [1] of an absorption peak and the transport mobility [44], respectively.

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.

2. 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]

$$\widehat{\mathcal{H}}{\mathsf{\Phi }}_{\mathbf{k}}\left(\mathit{r}\right)=E_{\mathbf{k}}{\mathsf{\Phi }}_{\mathbf{k}}\left(\mathit{r}\right),$$

(1)

where ${\mathsf{\Phi }}_{\mathbf{k}}\left(\mathit{r}\right)$ is the Bloch wave function, $E_{\mathbf{k}}$ the eigen-energy, and k is the wave vector of electrons within the first Brillouin zone of two-dimensional materials. The Hamiltonian operator $\widehat{\mathcal{H}}$ in Equation (1) takes a general form

$$\widehat{\mathcal{H}}={\widehat{\mathcal{H}}}_0+V\left(\mathit{r}\right),$$

(2)

in which the kinetic-energy operator ${\widehat{\mathcal{H}}}_0$ is

$${\widehat{\mathcal{H}}}_0=-\frac{{\hslash }^2}{2m_{\mathrm{e}}}{\nabla }_{\mathbf{r}}^2$$

(3)

with free-electron mass $m_{\mathrm{e}}$, while the potential energy $V\left(\mathit{r}\right)$ for an electron within the lattice of two-dimensional materials is given by [11]

$$V\left(\mathit{r}\right)=V_{\mathrm{ion}}\left(\mathit{r}\right)+\mathsf{\Delta }{V}_{\mathrm{L}}\left(\mathit{r}\right)$$

(4)

with $V_{\mathrm{ion}}\left(\mathit{r}\right)$ and $\mathsf{\Delta }{V}_{\mathrm{L}}\left(\mathit{r}\right)$ specifying the potentials of a single ion and that for the rest of ions, respectively. The Bloch wave function ${\mathsf{\Phi }}_{\mathbf{k}}\left(\mathit{r}\right)$ of electrons in Equation (1) can be decomposed into a linear combination of a set of orbital wave functions $\left\{{\varphi }_{\beta ,\mathbf{k}}\left(\mathit{r}\right)\right\}$ within the first Brillouin zone, leading to [11]

$${\mathsf{\Phi }}_{\mathbf{k}}\left(\mathit{r}\right)=\sum _{\beta }{C}_{\beta }{\varphi }_{\beta ,\mathbf{k}}\left(\mathit{r}\right),$$

(5)

where the index $\beta$ labels all the atomic orbitals of the lattice of two-dimensional materials. The expansion coefficient $C_{\beta }$ introduced in Equation (5) can be decided from

$$C_{\beta }=\int d^3\mathit{r}{\varphi }_{\beta ,\mathbf{k}}^{\ast }\left(\mathit{r}\right){\mathsf{\Phi }}_{\mathbf{k}}\left(\mathit{r}\right),$$

(6)

where the orthonormal property for the set of orbital wave functions $\left\{{\varphi }_{\beta ,\mathbf{k}}\left(\mathit{r}\right)\right\}$ 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 ${\varphi }_{\beta ,\mathbf{k}}\left(\mathit{r}\right)$ in Equation (5) by a linear combination of bond-orbital states $\left\{{\psi }_{\beta }(\mathit{r}-{\mathit{R}}_j)\right\}$ within a unit cell in real space, namely

$${\varphi }_{\beta ,\mathbf{k}}\left(\mathit{r}\right)=\frac{1}{\sqrt{N}}\sum _{j=1}^{N}exp(i\mathit{k}·{\mathit{R}}_j){\psi }_{\beta }(\mathit{r}-{\mathit{R}}_j),$$

(7)

where j is the index for all bonded lattice ions, ${\mathit{R}}_j$ the lattice-ion position vector, and N the total number of atoms within the unit cell. Here,

$${\psi }_{\beta }(\mathit{r}-{\mathit{R}}_j)=\frac{1}{\sqrt{N}}\sum _{\mathbf{k}}exp(-i\mathit{k}·{\mathit{R}}_j){\varphi }_{\beta ,\mathbf{k}}\left(\mathit{r}\right)$$

(8)

is termed as the localized Wannier function for the $\beta$ orbital of a bonded lattice ion at the site ${\mathit{R}}_j$, which satisfies the single-ion Schrödinger equation [16]

$$\left[{\widehat{\mathcal{H}}}_0+V_{\mathrm{ion}}\left(\mathit{r}\right)\right]{\psi }_{\beta }(\mathit{r}-{\mathit{R}}_j)={\epsilon }_{j,\beta }{\psi }_{\beta }(\mathit{r}-{\mathit{R}}_j)$$

(9)

with ${\epsilon }_{j,\beta }$ being the $\beta$th energy levels of electrons within an ion at the lattice site ${\mathit{R}}_j$.

Combining results in Equations (5) and (7), we acquire the following full LCAO expansion of a Bloch wave function [11]

$${\mathsf{\Phi }}_{\mathbf{k}}\left(\mathit{r}\right)=\frac{1}{\sqrt{N}}\sum _{j=1}^N\sum _{\beta }{C}_{\beta ;j}\left(\mathit{k}\right){\psi }_{\beta }(\mathit{r}-{\mathit{R}}_j)$$

(10)

with $C_{\beta ;j}\left(\mathit{k}\right)=C_{\beta }exp(i\mathit{k}·{\mathit{R}}_j)$. At the same time, using Equation (5), we find from Equation (1) that

$$\sum _{\beta }{C}_{\beta }\int d^2\mathit{r}{\varphi }_{\alpha ,\mathbf{k}}^{\ast }\left(\mathit{r}\right)\widehat{\mathcal{H}}{\varphi }_{\beta ,\mathbf{k}}\left(\mathit{r}\right)=E\left(\mathit{k}\right)\sum _{\beta }{C}_{\beta }\int d^2\mathit{r}{\varphi }_{\alpha ,\mathbf{k}}^{\ast }\left(\mathit{r}\right){\varphi }_{\beta ,\mathbf{k}}\left(\mathit{r}\right),$$

(11)

or equivalently, the following eigenvalue equation

$$\sum _{\beta }{\mathcal{H}}_{\alpha ,\beta }\left(\mathit{k}\right)C_{\beta }=E\left(\mathit{k}\right)\sum _{\beta }{\delta }_{\alpha ,\beta }{C}_{\beta }=E\left(\mathit{k}\right)C_{\alpha }.$$

(12)

As a result, the eigenvalue $E_n\left(\mathit{k}\right)$ can be determined from the secular determinant of Equation (12) for any given k, yielding

$$\mathcal{D}et\left\{{\mathcal{H}}_{\alpha ,\beta }\left(\mathit{k}\right)-E_n\left(\mathit{k}\right){\delta }_{\alpha ,\beta }\right\}=0,$$

(13)

and the orthonormal-eigenvectors $\left\{C_{n,\beta }\right\}$ are also obtained, corresponding to the eigenvalue $E_n\left(\mathit{k}\right)$ 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]

$${\mathcal{H}}_{\alpha \beta }\left(\mathit{k}\right)=\frac{1}{N}\sum _{i,j=1}^{N}exp\left(i\mathit{k}·{\mathit{R}}_j\right){\mathcal{H}}_{i\alpha ,j\beta },$$

(14)

in which ${\mathit{R}}_{ij}={\mathit{R}}_j-{\mathit{R}}_i$, and

$${\mathcal{H}}_{i\alpha ,j\beta }=\int d^2\mathit{r}{\psi }_{\alpha }^{\ast }(\mathit{r}-{\mathit{R}}_i)\left[{\widehat{\mathcal{H}}}_0+V_{\mathrm{ion}}\left(\mathit{r}\right)\right]{\psi }_{\beta }(\mathit{r}-{\mathit{R}}_j)+\int d^2\mathit{r}{\psi }_{\alpha }^{\ast }(\mathit{r}-{\mathit{R}}_i)\mathsf{\Delta }{V}_{\mathrm{L}}\left(\mathit{r}\right){\psi }_{\beta }(\mathit{r}-{\mathit{R}}_j).$$

(15)

In fact, we know from Equation (9) that

$$\int d^2\mathit{r}{\psi }_{\alpha }^{\ast }(\mathit{r}-{\mathit{R}}_i)\left[{\widehat{\mathcal{H}}}_0+V_{\mathrm{ion}}\left(\mathit{r}\right)\right]{\psi }_{\beta }(\mathit{r}-{\mathit{R}}_j)={\epsilon }_{j,\beta }{\delta }_{\alpha ,\beta }{\delta }_{i,j},$$

(16)

$$\int d^2\mathit{r}{\psi }_{\alpha }^{\ast }(\mathit{r}-{\mathit{R}}_i)\mathsf{\Delta }{V}_{\mathrm{L}}\left(\mathit{r}\right){\psi }_{\beta }(\mathit{r}-{\mathit{R}}_j)\equiv \left\{\begin{array}{cc}{C}_{\mathsf{\Sigma }}{\delta }_{\alpha ,\beta }& \mathrm{if} i=j,\\ \\ t_{\alpha \beta }\left({\mathit{R}}_{ij}\right)& \mathrm{if} i\ne j,\end{array}\right.$$

(17)

where $({\epsilon }_{j,\beta }+C_{\mathsf{\Sigma }})$ represents the site energy, and $t_{\alpha \beta }\left({\mathit{R}}_{ij}\right)$ 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

$$\begin{array}{l}\int d^2\mathit{r}{\mathsf{\Phi }}_{n^{\prime },{\mathbf{k}}^{\prime }}^{\ast }\left(\mathit{r}\right)\widehat{\mathcal{H}}{\mathsf{\Phi }}_{n,\mathbf{k}}\left(\mathit{r}\right)=\frac{1}{N}{\displaystyle \sum _{j,j^{\prime }=1}^N}{\displaystyle \sum _{\alpha ,\beta }}{C}_{n^{\prime },\alpha }^{\ast }{C}_{n,\beta }exp(i\mathit{k}·{\mathit{R}}_j-i{\mathit{k}}^{\prime }·{\mathit{R}}_{j^{\prime }}){\mathcal{H}}_{j^{\prime }\alpha ,j\beta }\hfill \\ =\frac{({\epsilon }_{j,\beta }+C_{\mathsf{\Sigma }})}{N}{\displaystyle \sum _{j=1}^N}{\displaystyle \sum _{\beta }}{C}_{n^{\prime },\beta }^{\ast }{C}_{n,\beta }exp[i(\mathit{k}-{\mathit{k}}^{\prime })·{\mathit{R}}_j]\hfill \\ +\frac{1}{N}{\displaystyle \sum _{j,j^{\prime }=1}^N}{}^{\prime }{\displaystyle \sum _{\alpha ,\beta }}{C}_{n^{\prime },\alpha }^{\ast }{C}_{n,\beta }exp(i\mathit{k}·{\mathit{R}}_j-i{\mathit{k}}^{\prime }·{\mathit{R}}_{j^{\prime }})t_{\alpha \beta }\left({\mathit{R}}_{j^{\prime }j}\right),\hfill \end{array}$$

(18)

where the primed summation in the second term of the right-hand side of the last equation excludes the contribution from $j=j^{\prime }$, and $C_{n,\beta }$ can be obtained from the calculated eigenvector from Equation (12). The matrix elements for other physical operators can be computed in a similar way.

3. 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,\cdots$ will degenerate with the angular-momentum quantum number $\ell =0,1,\cdots ,n-1$, as well as the magnetic quantum number $m=-\ell ,\cdots ,0,\cdots ,\ell$ [52]. Consequently, there exists a total orbital degeneracy $n^2$ (excluding the spin-degeneracy). Customarily, we specify these orbitals by $\ell =0,1,2,3,\cdots$ for $\{s,p,d,f,\cdots \}$ 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,\cdots \}$. To further specify the spatial direction of the chemical bonding between two atoms at the lattice sites ${\mathit{R}}_i$ and ${\mathit{R}}_j$, we have to rely on three directional cosines $\ell ,m,n$, as defined in Figure 2.

Considering s and p orbitals as an example, we display their possible bonding potentials $V_{\ell ,{\ell }^{\prime };\sigma \left(\pi \right)}$ in Figure 3 for s, p orbitals and four different configurations, including $\sigma$ and $\pi$ bonds. Meanwhile, we also list six different $\pi$, $\sigma$, $\delta$ bonding configurations in Figure 4 for s, p, d orbitals.

To speed up numerical computations, the bonding potentials $V_{\ell ,{\ell }^{\prime };\sigma \left(\pi \right)}$ for $\ell ,{\ell }^{\prime }=s,p,d$ in Figure 3 and Figure 4 are usually parameterized as: [53] $V_{\ell ,{\ell }^{\prime };\xi }=({\hslash }^2/m_{\mathrm{e}}{d}^2){\eta }_{\ell {\ell }^{\prime };\xi }$, $V_{\ell ,d;\xi }=({\hslash }^2{r}_d^{3/2}/m_{\mathrm{e}}{d}^{7/2}){\eta }_{\ell ,d;\xi }$ and $V_{d,d;\xi }=({\hslash }^2{r}_d^3/m_{\mathrm{e}}{d}^5){\eta }_{d,d;\xi }$, where d and $r_d$ represent the bonding length and atomic radius, and $\xi =\sigma ,\pi ,\delta$ are for various bond configurations. Here, the dimensionless bonding parameters ${\eta }_{\ell ,{\ell }^{\prime };\xi }$ for different bonding types are listed in

Table 1. Inter-atomic bonding parameters.

Bonding Parameter	Value [11]
${\eta }_{s,s;\sigma }$	$-1.40$
${\eta }_{s,d;\sigma }$	$-3.16$
${\eta }_{d,d;\sigma }$	$-16.2$
${\eta }_{s,p;\sigma }$	$1.84$
${\eta }_{p,d;\sigma }$	$-2.95$
${\eta }_{d,d;\pi }$	$8.75$
${\eta }_{p,p;\sigma }$	$3.24$
${\eta }_{p,d;\pi }$	$1.36$
${\eta }_{d,d;\delta }$	0
${\eta }_{p,p;\pi }$	$-0.81$

.

By using these parameterized bonding potentials $V_{\ell ,{\ell }^{\prime };\xi }$, $V_{\ell ,d;\xi }$ and $V_{d,d;\xi }$, we are able to compute further the hopping integrals $t_{\alpha \beta }\left({\mathit{R}}_{ij}\right)$ based on the Slater-Koster approximation [14], and some commonly-used results are shown in

Table 2. Expressions for Bonding Integrals.

Bonding Integral	Expression [14]
$t_{s,s}$	$V_{s,s;\sigma }$
$t_{s,x}$	$\ell V_{s,p;\sigma }$
$t_{x,x}$	${\ell }^2{V}_{p,p;\sigma }+(1-{\ell }^2)V_{p,p;\pi }$
$t_{x,y}$	$\ell m(V_{p,p;\sigma }-V_{p,p;\pi })$
$t_{x,z}$	$\ell n(V_{p,p;\sigma }-V_{p,p;\pi })$
$t_{s,xy}$	$\sqrt{3}\ell m{V}_{s,d;\sigma }$
$t_{s,x^2-y^2}$	$(\sqrt{3}/2)({\ell }^2-m^2)V_{s,d;\sigma }$
$t_{3z^2-r^2}$	$[(n^2-({\ell }^2+m^2)/2]V_{s,d;\sigma }$
⋮	⋮

.

4. Tight-Binding Model for Graphene Band Structure

To seek for an application, we use the general theory, as developed in Section 2 and Section 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 $\sigma$ 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 $\sigma$ bonds, while the remaining electron in the z orbital leads to the lower-energy $\pi$ bonds, i.e., a side-on overlap of the $2p$-orbital wave functions. Consequently, these $\pi$-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 $\pi$-bond ($p_z$-orbital) electrons in graphene can be expressed as

$${\varphi }_{p_z,\mathbf{k}}^{A,B}\left(\mathit{r}\right)=\frac{1}{\sqrt{N}}\sum _{j=1}^{N}exp(i\mathit{k}·{\mathit{R}}_j){\psi }_{p_z}^{A,B}(\mathit{r}-{\mathit{R}}_j),$$

(19)

where $\mathit{k}\equiv (k_x,k_y)$, ${\mathit{R}}_j\equiv (R_j^x,R_j^y)=m_j{\mathit{a}}_1+n_j{\mathit{a}}_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

$${\varphi }_{p_z,\mathbf{k}}\left(\mathit{r}\right)=a_{\mathbf{k}}{\varphi }_{p_z,\mathbf{k}}^A\left(\mathit{r}\right)+b_{\mathbf{k}}{\varphi }_{p_z,\mathbf{k}}^B\left(\mathit{r}\right),$$

(20)

where $a_{\mathbf{k}}$ and $b_{\mathbf{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

$$\left[\begin{array}{cc}{\mathcal{H}}_{AA}\left(\mathit{k}\right)& {\mathcal{H}}_{AB}\left(\mathit{k}\right)\\ \\ {\mathcal{H}}_{BA}\left(\mathit{k}\right)& {\mathcal{H}}_{BB}\left(\mathit{k}\right)\end{array}\right]\left[\begin{array}{c}{a}_{\mathbf{k}}\\ \\ b_{\mathbf{k}}\end{array}\right]=E_n\left(\mathit{k}\right)\left[\begin{array}{cc}{\mathcal{S}}_{AA}\left(\mathit{k}\right)& {\mathcal{S}}_{AB}\left(\mathit{k}\right)\\ \\ {\mathcal{S}}_{BA}\left(\mathit{k}\right)& {\mathcal{S}}_{BB}\left(\mathit{k}\right)\end{array}\right]\left[\begin{array}{c}{a}_{\mathbf{k}}\\ \\ b_{\mathbf{k}}\end{array}\right],$$

(21)

where ${\mathcal{S}}_{\ell {\ell }^{\prime }}\left(\mathit{k}\right)=\langle {\varphi }_{p_z,\mathbf{k}}^{\ell }|{\varphi }_{p_z,\mathbf{k}}^{{\ell }^{\prime }}\rangle$, ${\mathcal{H}}_{\ell {\ell }^{\prime }}\left(\mathit{k}\right)=\langle {\varphi }_{p_z,\mathbf{k}}^{\ell }|\widehat{\mathcal{H}}|{\varphi }_{p_z,\mathbf{k}}^{{\ell }^{\prime }}\rangle$, $\ell ,{\ell }^{\prime }=$ A or B, and $E_n\left(\mathit{k}\right)$ represents the eigen-energies of $\pi$-bond electrons with $n=1,2$ labeling two graphene low-energy bands determined by the secular determinant: $\mathcal{D}et{\left\{{\mathcal{H}}_{\ell {\ell }^{\prime }}\left(\mathit{k}\right)-E_n\left(\mathit{k}\right){\mathcal{S}}_{\ell {\ell }^{\prime }}\left(\mathit{k}\right)\right\}}_{2\times 2}=0$.

As in Equation (7), we can rewrite the orbital wave function ${\varphi }_{p_z,\mathbf{k}}\left(\mathit{r}\right)$ in Equation (20) approximately only by its near-neighbor decomposition, yielding

$${\varphi }_{p_z,\mathbf{k}}\left(\mathit{r}\right)=\frac{1}{\sqrt{N_c}}\sum _{\ell \in A,B}exp(i\mathit{k}·{\mathit{R}}_{\ell })\sum _{j=1}^{N_c}{a}_{j\mathbf{k}}{\psi }_{p_z}(\mathit{r}-{\mathit{R}}_{\ell }+{\mathsf{\Delta }}_j),$$

(22)

and then, the eigenvalue equation turns into $\mathcal{D}et{\left\{{\mathcal{H}}_{j{j}^{\prime }}\left(\mathit{k}\right)-E_s\left(\mathit{k}\right){\mathcal{S}}_{j{j}^{\prime }}\left(\mathit{k}\right)\right\}}_{N_c\times N_c}=0$ with eigen-vectors ${\left\{a_{j\mathbf{k}}\right\}}_{N_c\times 1}$, where $N_c$ represents the number of near-neighbor atoms within a unit cell, ${\mathsf{\Delta }}_j$ stands for the lattice vectors of the near-neighbor atoms relative to the sublattice site ${\mathit{R}}_{\ell }$, and $a_{j\mathbf{k}}=a_{\mathbf{k}}exp(-i\mathit{k}·{\mathsf{\Delta }}_j)$. Moreover, we find ${\mathcal{H}}_{j{j}^{\prime }}\left(\mathit{k}\right)={\epsilon }_2{\mathcal{S}}_{j{j}^{\prime }}\left(\mathit{k}\right)+t_{j{j}^{\prime }}\left(\mathit{k}\right)$, where ${\epsilon }_2$ stands for the second energy level of electrons within a carbon atom,

$${\mathcal{S}}_{j{j}^{\prime }}\left(\mathit{k}\right)=\sum _{\ell \in A,B}exp(i\mathit{k}·{\mathit{R}}_{\ell })\int d^2\mathit{r}{\psi }_{p_z}^{\ast }(\mathit{r}-{\mathit{R}}_{\ell }+{\mathsf{\Delta }}_j){\psi }_{p_z}(\mathit{r}-{\mathit{R}}_{\ell }+{\mathsf{\Delta }}_{j^{\prime }})$$

(23)

is the overlap integral, while

$$t_{j{j}^{\prime }}\left(\mathit{k}\right)=\sum _{\ell \in A,B}exp(i\mathit{k}·{\mathit{R}}_{\ell })\int d^2\mathit{r}{\psi }_{p_z}^{\ast }(\mathit{r}-{\mathit{R}}_{\ell }+{\mathsf{\Delta }}_j)\mathsf{\Delta }{V}_L\left(\mathit{r}\right){\psi }_{p_z}(\mathit{r}-{\mathit{R}}_{\ell }+{\mathsf{\Delta }}_{j^{\prime }})$$

(24)

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 ${\mathbf{\delta }}_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

$$\begin{array}{ccc}\hfill t_{AB}\left(\mathit{k}\right)& =& t_{BA}^{\ast }\left(\mathit{k}\right)={\gamma }^{\ast }\left(\mathit{k}\right)t_{nn},\hfill \\ \hfill t_{AA}\left(\mathit{k}\right)-C_p{\mathcal{S}}_{AA}& =& t_{BB}\left(\mathit{k}\right)-C_p{\mathcal{S}}_{BB}=2t_{nnn}\sum _{i=1}^{3}cos(\mathit{k}·{\mathit{a}}_i)=\left(\right|\gamma \left(\mathit{k}\right){|}^2-3)t_{nnn},\hfill \\ \hfill {\mathcal{S}}_{AB}\left(\mathit{k}\right)& =& S_{BA}^{\ast }\left(\mathit{k}\right)={\gamma }^{\ast }\left(\mathit{k}\right)s_{nn},\hfill \\ \hfill {\mathcal{S}}_{AA}\left(\mathit{k}\right)& =& {\mathcal{S}}_{BB}\left(\mathit{k}\right)=1+{\left(\right|\gamma \left(\mathit{k}\right)|}^2-3)s_{nnn}\approx 1,\hfill \end{array}$$

(25)

where ${\mathit{a}}_3\equiv {\mathit{a}}_1-{\mathit{a}}_2$, $\gamma \left(\mathit{k}\right)=1+exp(i\mathit{k}·{\mathit{a}}_1)+exp(i\mathit{k}·{\mathit{a}}_2)$, and the hopping and overlap integrals are calculated as

$$\begin{array}{ccc}\hfill C_p& =& \int d^2\mathit{r}{\psi }_A^{\ast }\left(\mathit{r}\right)\mathsf{\Delta }{V}_L\left(\mathit{r}\right){\psi }_A\left(\mathit{r}\right)=\int d^2\mathit{r}{\psi }_B^{\ast }\left(\mathit{r}\right)\mathsf{\Delta }{V}_L\left(\mathit{r}\right){\psi }_B\left(\mathit{r}\right),\hfill \\ \hfill t_{nn}& =& \int d^2\mathit{r}{\psi }_A^{\ast }\left(\mathit{r}\right)\mathsf{\Delta }{V}_L\left(\mathit{r}\right){\psi }_B(\mathit{r}+{\mathbf{\delta }}_3),\hfill \\ \hfill t_{nnn}& =& \int d^2\mathit{r}{\psi }_A^{\ast }\left(\mathit{r}\right)\mathsf{\Delta }{V}_L\left(\mathit{r}\right){\psi }_A(\mathit{r}+{\mathit{a}}_1)=\int d^2\mathit{r}{\psi }_B^{\ast }\left(\mathit{r}\right)\mathsf{\Delta }{V}_L\left(\mathit{r}\right){\psi }_B(\mathit{r}+{\mathit{a}}_1),\hfill \\ \hfill s_{nn}& =& \int d^2\mathit{r}{\psi }_A^{\ast }\left(\mathit{r}\right){\psi }_B(\mathit{r}+{\mathbf{\delta }}_3),\hfill \\ \hfill s_{nnn}& =& \int d^2\mathit{r}{\psi }_A^{\ast }\left(\mathit{r}\right){\psi }_A(\mathit{r}+{\mathit{a}}_1)=\int d^2\mathit{r}{\psi }_B^{\ast }\left(\mathit{r}\right){\psi }_B(\mathit{r}+{\mathit{a}}_1).\hfill \end{array}$$

(26)

Particularly, the results for these tight-binding model parameters in Equation (26) for band structures are presented in

Table 3. Graphene structure and tight-binding model parameters.

Parameter	Value [55]
${\mathit{a}}_1$	$(a/2)(3,\sqrt{3})$
${\mathit{a}}_2$	$(a/2)(3,-\sqrt{3})$
${\mathit{a}}_3$	$(a/2)(0,2\sqrt{3})$
${\mathbf{\delta }}_1$	$(a/2)(1,\sqrt{3})$
${\mathbf{\delta }}_2$	$(a/2)(1,-\sqrt{3})$
${\mathbf{\delta }}_3$	$-a(1,0)$
K	$(2\pi /3\sqrt{3}a)(\sqrt{3},1)$
${\mathit{K}}^{\prime }$	$(2\pi /3\sqrt{3}a)(\sqrt{3},-1)$
$s_{nn}$	$0.106$
$s_{nnn}$	$0.001$
$t_{nn}$	$-2.78$eV
$t_{nnn}$	$-0.12$eV
${\epsilon }_2+C_p$	$-0.36$eV

, which have been computed from listed bonding parameters in Table 1 and bonding integrals in Table 2.

Finally, from the eigenvalue equation $\mathcal{D}et{\left\{t_{\ell {\ell }^{\prime }}\left(\mathit{k}\right)-E\left(\mathit{k}\right){\mathcal{S}}_{\ell {\ell }^{\prime }}\left(\mathit{k}\right)\right\}}_{2\times 2}=0$ in Equation (21) for $\ell ,{\ell }^{\prime }=A,B$, we obtain an explicit expression

$$E^2\left(\mathit{k}\right)\mathcal{D}et\left\{\overleftrightarrow{\mathcal{S}}\right\}-E\left(\mathit{k}\right)\left[{\mathcal{S}}_{AA}{t}_{BB}+{\mathcal{S}}_{BB}{t}_{AA}-{\mathcal{S}}_{AB}{t}_{BA}-{\mathcal{S}}_{BA}{t}_{AB}\right]+\mathcal{D}et\left\{t\right\}=0,$$

(27)

where, by setting $s_{nnn}=0$, we have three coefficients

$$\begin{array}{ccc}\hfill & & \mathcal{D}et\left\{\overleftrightarrow{\mathcal{S}}\right\}=1-s_{nn}^2{\left|\gamma \left(\mathit{k}\right)\right|}^2,\hfill \\ \hfill & & \mathcal{D}et\left\{\overleftrightarrow{t}\right\}=(|\gamma \left(\mathit{k}\right){|}^2-{3\left|\right)}^2{t}_{nnn}^2-t_{nn}^2{\left|\gamma \left(\mathit{k}\right)\right|}^2,\hfill \\ & & {\mathcal{S}}_{AA}{t}_{BB}+{\mathcal{S}}_{BB}{t}_{AA}-{\mathcal{S}}_{AB}{t}_{BA}-{\mathcal{S}}_{BA}{t}_{AB}=2\left[{\left(\right|\gamma \left(\mathit{k}\right)|}^2-3)t_{nnn}-t_{nn}{s}_{nn}{\left|\gamma \left(\mathit{k}\right)\right|}^2\right].\hfill \end{array}$$

(28)

This leads to the explicit solution of Equation (27), namely [55]

$$E_{\lambda }\left(\mathit{k}\right)=\left({\epsilon }_2+C_p\right)+\frac{{\left(\right|\gamma \left(\mathit{k}\right)|}^2-3\left|\right)t_{nnn}-t_{nn}{s}_{nn}{\left|\gamma \left(\mathit{k}\right)\right|}^2+\lambda \sqrt{\mathcal{D}\left(\mathit{k}\right)}}{1-s_{nn}^2{\left|\gamma \left(\mathit{k}\right)\right|}^2},$$

(29)

where $\lambda =\pm 1$ correspond to valence ($-1$) and conduction ($+1$) bands, respectively, and

$$\begin{array}{ccc}\hfill \mathcal{D}\left(\mathit{k}\right)& =& {\left[{\left(\right|\gamma \left(\mathit{k}\right)|}^2-3)t_{nnn}-t_{nn}{s}_{nn}{\left|\gamma \left(\mathit{k}\right)\right|}^2\right]}^2-(1-s_{nn}^2{\left|\gamma \left(k\right)\right|}^2)\hfill \\ & \times & \left[{\left(\right|\gamma \left(\mathit{k}\right)|}^2-{3\left|\right)}^2{t}_{nnn}^2-t_{nn}^2{\left|\gamma \left(\mathit{k}\right)\right|}^2\right]={\left|\gamma \left(k\right)\right|}^2{\left[{\left(\right|\gamma \left(\mathit{k}\right)|}^2-3)t_{nnn}{s}_{nn}+t_{nn}\right]}^2.\hfill \end{array}$$

(30)

By using the result in Equations (29) and (30) can be rewritten as

$$\begin{array}{ccc}\hfill E_{\lambda }\left(\mathit{k}\right)& =& \left({\epsilon }_2+C_p\right)+\frac{{\left(\right|\gamma \left(\mathit{k}\right)|}^2-3)t_{nnn}[1+\lambda |\gamma \left(\mathit{k}\right)|s_{nn}]-t_{nn}{s}_{nn}{\left|\gamma \left(\mathit{k}\right)\right|}^2+\lambda \left|\gamma \left(\mathit{k}\right)\right|t_{nn}}{1-s_{nn}^2{\left|\gamma \left(\mathit{k}\right)\right|}^2}\hfill \\ & \approx & \left({\epsilon }_2+C_p\right)+\frac{{\left(\right|\gamma \left(\mathit{k}\right)|}^2-3)t_{nnn}-t_{nn}{s}_{nn}{\left|\gamma \left(\mathit{k}\right)\right|}^2+\lambda t_{nn}\left|\gamma \left(\mathit{k}\right)\right|}{1-s_{nn}^2{\left|\gamma \left(\mathit{k}\right)\right|}^2}.\hfill \end{array}$$

(31)

By setting $C_p+{\epsilon }_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_{\mathbf{k}}^{\lambda }$ and $b_{\mathbf{k}}^{\lambda }$, are found to be

$$\begin{array}{ccc}\hfill a_{\mathbf{k}}^{\lambda }& =& \left\{\frac{{\gamma }^{\ast }\left(\mathit{k}\right)[E_{\lambda }\left(\mathit{k}\right)s_{nn}-t_{nn}]}{{\left(\right|\gamma \left(k\right)|}^2-3)t_{nnn}-E_{\lambda }\left(\mathit{k}\right)}\right\}{b}_{\mathbf{k}}^{\lambda },\hfill \\ \hfill b_{\mathbf{k}}^{\lambda }& =& \frac{{\left|\right(\left|\gamma \left(\mathit{k}\right)\right|}^2-3)t_{nnn}-E_{\lambda }\left(\mathit{k}\right)|}{\sqrt{|{\gamma }^{\ast }\left(\mathit{k}\right)[t_{nn}-E_{\lambda }\left(\mathit{k}\right)s_{nn}]{|}^2+{\left|\right(\left|\gamma \left(k\right)\right|}^2-3{)t_{nnn}-E_{\lambda }\left(k\right)|}^2}}.\hfill \end{array}$$

(32)

As known experimentally, both the nearest-neighbor (nn) overlap and the next-nearest-neighbor (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

$$E_{\lambda }\left(\mathit{k}\right)\approx 2t_{nnn}^{\prime }\sum _{i=1}^{3}cos(\mathit{k}·{\mathit{a}}_i)+\lambda t_{nn}{\left[3+2\sum _{i=1}^{3}cos(\mathit{k}·{\mathit{a}}_i)\right]}^{1/2},$$

(33)

where $t_{nnn}^{\prime }=t_{nnn}-s_{nn}{t}_{nn}$ is the corrected hopping amplitude.

5. 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 tight-binding 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_{\mathrm{e}}=n_{\mathrm{h}}\equiv 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 ${\mathsf{\Delta }}_{\mathrm{e}}\left(\mathit{k}\right)$ and ${\mathsf{\Delta }}_{\mathrm{h}}\left(\mathit{k}\right)$ for electrons (e) and holes (h), respectively. These two rates determine the inhomogeneous line-shape of a resonant interband-absorption peak at $\hslash \omega ={\epsilon }_{\mathbf{k}}^{\mathrm{e}}+{\epsilon }_{\mathbf{k}}^{\mathrm{h}}$ for vertical transitions of electrons with their kinetic energies ${\epsilon }_{\mathbf{k}}^{\mathrm{e},\mathrm{h}}$ in valence and conduction bands.

As illustrated by Feynman diagrams [3] in Figure 7, the CDD rate ${\mathsf{\Delta }}_{\mathrm{e}}\left(\mathit{k}\right)$ of electrons is calculated as [1,44]

$$\begin{array}{ccc}\hfill {\mathsf{\Delta }}_{\mathrm{e}}\left(\mathit{k}\right)& =& \frac{8\pi }{\hslash {\mathcal{A}}^2}\sum _{{\mathbf{k}}_1,\mathbf{q}\ne 0}{\left|V_{\mathbf{k},{\mathbf{k}}_1;{\mathbf{k}}_1-\mathbf{q},\mathbf{k}+\mathbf{q}}^{\mathrm{ee}}\right|}^2\left[{\mathcal{L}}_0({\epsilon }_{{\mathbf{k}}_1-\mathbf{q}}^{\mathrm{e}}+{\epsilon }_{\mathbf{k}+\mathbf{q}}^{\mathrm{e}}-{\epsilon }_{{\mathbf{k}}_1}^{\mathrm{e}}-{\epsilon }_{\mathbf{k}}^{\mathrm{e}},{\mathsf{\Gamma }}_{\mathrm{e}})\right.\hfill \\ \hfill & \times & \left.\left\{f_{{\mathbf{k}}_1-\mathbf{q}}^{\mathrm{e}}{f}_{\mathbf{k}+\mathbf{q}}^{\mathrm{e}}(1-f_{{\mathbf{k}}_1}^{\mathrm{e}})+(1-f_{{\mathbf{k}}_1-\mathbf{q}}^{\mathrm{e}})(1-f_{\mathbf{k}+\mathbf{q}}^{\mathrm{e}})f_{{\mathbf{k}}_1}^{\mathrm{e}}\right\}\right]\hfill \\ \hfill & +& \frac{8\pi }{\hslash {\mathcal{A}}^2}\sum _{{\mathbf{k}}_1,\mathbf{q}\ne 0}{\left|V_{\mathbf{k},{\mathbf{k}}_1;{\mathbf{k}}_1-\mathbf{q},\mathbf{k}-\mathbf{q}}^{\mathrm{he}}\right|}^2\left[{\mathcal{L}}_0({\epsilon }_{{\mathbf{k}}_1-\mathbf{q}}^{\mathrm{e}}+{\epsilon }_{-(\mathbf{k}-\mathbf{q})}^{\mathrm{h}}-{\epsilon }_{{\mathbf{k}}_1}^{\mathrm{e}}-{\epsilon }_{-\mathbf{k}}^{\mathrm{h}},{\mathsf{\Gamma }}_{\mathrm{eh}})\right.\hfill \\ & \times & \left.\left\{f_{{\mathbf{k}}_1-\mathbf{q}}^{\mathrm{e}}{f}_{-(\mathbf{k}-\mathbf{q})}^{\mathrm{h}}(1-f_{{\mathbf{k}}_1}^{\mathrm{e}})+(1-f_{{\mathbf{k}}_1-\mathbf{q}}^{\mathrm{e}})(1-f_{-(\mathbf{k}-\mathbf{q})}^{\mathrm{h}})f_{{\mathbf{k}}_1}^{\mathrm{e}}\right\}\right],\hfill \end{array}$$

(34)

where both spin and valley degeneracies are included, $\mathcal{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, ${\epsilon }_{\mathbf{k}}^{\mathrm{e},\mathrm{h}}$ in Equation (34) stands for the kinetic energy of electrons (e) or holes (h), and $f_{\mathbf{k}}^{\mathrm{e},\mathrm{h}}={\{1+exp[({\epsilon }_{\mathbf{k}}^{\mathrm{e},\mathrm{h}}-{\mu }_{\mathrm{e},\mathrm{h}})/k_{\mathrm{B}}T]\}}^{-1}$ is the Fermi function for thermal-equilibrium photo-carriers with their chemical potentials ${\mu }_{\mathrm{e},\mathrm{h}}$ at temperature T. Here, ${\mu }_{\mathrm{e},\mathrm{h}}$ are separately determined by following two equations for given T, i.e.,

$$n_{\mathrm{e},\mathrm{h}}=\frac{4}{\mathcal{A}}\sum _{\mathbf{k}}\frac{1}{1+exp[({\epsilon }_{\mathbf{k}}^{\mathrm{e},\mathrm{h}}-{\mu }_{\mathrm{e},\mathrm{h}})/k_{\mathrm{B}}T]},$$

(35)

where both spin and valley degeneracies are included and ${\mu }_{\mathrm{e}}={\mu }_{\mathrm{h}}$ in our case. Furthermore, in Equation (34), ${\mathcal{L}}_0(a,b)=(b/\pi )/(a^2+b^2)$ is the Lorentzian line-shape function, ${\mathsf{\Gamma }}_{\mathrm{e},\mathrm{h}}$ are inverse lifetime of unperturbed electrons or holes, and ${\mathsf{\Gamma }}_{\mathrm{eh}}=({\mathsf{\Gamma }}_{\mathrm{e}}+{\mathsf{\Gamma }}_{\mathrm{h}})/2$.

In addition, we have introduced in Equation (34), as well as in Equation (40) below, the Coulomb-interaction matrix elements, given by [56]

$$\begin{array}{ccc}\hfill V_{\mathbf{k},{\mathbf{k}}_1;{\mathbf{k}}_1-\mathbf{q},\mathbf{k}+\mathbf{q}}^{\mathrm{ee}}& =& u_{\mathrm{c}}\left(q\right){\mathcal{F}}_{\mathbf{k},\mathbf{k}+\mathbf{q}}^{\left(\mathrm{c}\right)}\left(\mathit{q}\right){\mathcal{F}}_{{\mathbf{k}}_1,{\mathbf{k}}_1-\mathbf{q}}^{\left(\mathrm{c}\right)}(-\mathit{q}),\hfill \\ \hfill V_{\mathbf{k},{\mathbf{k}}_1;{\mathbf{k}}_1-\mathbf{q},\mathbf{k}-\mathbf{q}}^{\mathrm{he}}& =& u_{\mathrm{c}}\left(q\right){\mathcal{F}}_{\mathbf{k},\mathbf{k}-\mathbf{q}}^{\left(\mathrm{v}\right)}\left(\mathit{q}\right){\mathcal{F}}_{{\mathbf{k}}_1,{\mathbf{k}}_1-\mathbf{q}}^{\left(\mathrm{c}\right)}(-\mathit{q}),\hfill \\ \hfill V_{\mathbf{k},{\mathbf{k}}_1;{\mathbf{k}}_1-\mathbf{q},\mathbf{k}+\mathbf{q}}^{\mathrm{hh}}& =& u_{\mathrm{c}}\left(q\right){\mathcal{F}}_{\mathbf{k},\mathbf{k}+\mathbf{q}}^{\left(\mathrm{v}\right)}\left(\mathit{q}\right){\mathcal{F}}_{{\mathbf{k}}_1,{\mathbf{k}}_1-\mathbf{q}}^{\left(\mathrm{v}\right)}(-\mathit{q}),\hfill \\ \hfill V_{\mathbf{k},{\mathbf{k}}_1;{\mathbf{k}}_1-\mathbf{q},\mathbf{k}-\mathbf{q}}^{\mathrm{eh}}& =& u_{\mathrm{c}}\left(q\right){\mathcal{F}}_{\mathbf{k},\mathbf{k}-\mathbf{q}}^{\left(\mathrm{c}\right)}\left(\mathit{q}\right){\mathcal{F}}_{{\mathbf{k}}_1,{\mathbf{k}}_1-\mathbf{q}}^{\left(\mathrm{v}\right)}(-\mathit{q}),\hfill \end{array}$$

(36)

where $u_{\mathrm{c}}\left(q\right)=e^2/\left[2{ϵ}_0{ϵ}_{\mathrm{r}}(q+q_0)\right])$ in Equation (36) is the two-dimensional Fourier transformed Coulomb potential $\sim 1/r$ including static screening, ${ϵ}_0$ represents the vacuum permittivity, and ${ϵ}_{\mathrm{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]

$$q_0=\left(\frac{e^2}{8{ϵ}_0{ϵ}_{\mathrm{r}}{k}_{\mathrm{B}}T}\right)\frac{4}{\mathcal{A}}\sum _{\mathbf{k}}\left[{cosh}^{-2}\left(\frac{{\epsilon }_{\mathbf{k}}^{\mathrm{e}}-{\mu }_{\mathrm{e}}}{2k_{\mathrm{B}}T}\right)+{cosh}^{-2}\left(\frac{{\epsilon }_{\mathbf{k}}^{\mathrm{h}}-{\mu }_{\mathrm{h}}}{2k_{\mathrm{B}}T}\right)\right],$$

(37)

where both spin and valley degeneracies have been included.

Furthermore, the introduced ${\mathcal{F}}_{\mathbf{k},{\mathbf{k}}^{\prime }}^{\left(s\right)}\left(q\right)$ in Equation (36) with $s=c,v$ represents the Bloch-function form factor, calculated as [57]

$$\begin{array}{ccc}\hfill & & {\mathcal{F}}_{\mathbf{k},{\mathbf{k}}^{\prime }}^{\left(s\right)}\left(\mathit{q}\right)=\int d^2\mathit{r}{\left[{\mathsf{\Phi }}_{\mathbf{k}}^{\left(s\right)}\left(\mathit{r}\right)\right]}^{\ast }exp\left(i\mathit{q}·\mathit{r}\right){\mathsf{\Phi }}_{{\mathbf{k}}^{\prime }}^{\left(s\right)}\left(\mathit{r}\right)=\frac{1}{N_c}\left\{{\left[{\mathit{a}}_{\mathbf{k}}^{\left(s\right)}\right]}^{†}\otimes {\mathit{a}}_{{\mathbf{k}}^{\prime }}^{\left(s\right)}\right\}\hfill \\ & \times & \sum _{j,j^{\prime }=1}^{N_c}exp[-i(\mathit{q}-\mathit{k})·{\mathsf{\Delta }}_j-i{\mathit{k}}^{\prime }·{\mathsf{\Delta }}_{j^{\prime }}]\sum _{\ell ,{\ell }^{\prime }\in A,B}exp[i({\mathit{k}}^{\prime }-\mathit{k}+\mathit{q})·{\mathit{R}}_{\ell }]{\mathcal{W}}_s(\mathit{q},{\mathit{R}}_{{\ell }^{\prime }\ell }+{\mathsf{\Delta }}_{j{j}^{\prime }}),\hfill \end{array}$$

(38)

where the Bloch functions ${\mathsf{\Phi }}_{\mathbf{k}}^{c,v}\left(\mathit{r}\right)$ in Equations (5) and (22) have been employed. In Equation (38), $N_c$ represents the number of near-neighbor atoms within a unit cell, ${\mathsf{\Delta }}_j$ stands for the lattice vectors of the near-neighbor atoms relative to the sublattice site ${\mathit{R}}_{\ell }$, and ${\mathit{a}}_{\mathbf{k}}^{\left(s\right)}$ are two column eigenvectors in Equation (32) for $s=c,v$. The Wannier-function structure factor ${\mathcal{W}}_s(\mathit{q},{\mathit{R}}_{{\ell }^{\prime }\ell }+{\mathsf{\Delta }}_{j{j}^{\prime }})$ in Equation (38) is defined as

$${\mathcal{W}}_s(\mathit{q},{\mathit{R}}_{{\ell }^{\prime }\ell }+{\mathsf{\Delta }}_{j{j}^{\prime }})=\int d^2\mathit{r}{\left[{\psi }_{p_z}^{\left(s\right)}\left(\mathit{r}\right)\right]}^{\ast }exp(i\mathit{q}·\mathit{r}){\psi }_{p_z}^{\left(s\right)}(\mathit{r}-{\mathit{R}}_{{\ell }^{\prime }\ell }-{\mathsf{\Delta }}_{j{j}^{\prime }}),$$

(39)

where ${\mathit{R}}_{{\ell }^{\prime }\ell }={\mathit{R}}_{{\ell }^{\prime }}-{\mathit{R}}_{\ell }$ and ${\mathsf{\Delta }}_{j{j}^{\prime }}={\mathsf{\Delta }}_j-{\mathsf{\Delta }}_{j^{\prime }}$. 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 ${\mathsf{\Delta }}_{\mathrm{h}}\left(\mathit{k}\right)$ of holes takes the form [1,44]

$$\begin{array}{ccc}\hfill {\mathsf{\Delta }}_{\mathrm{h}}\left(\mathit{k}\right)& =& \frac{8\pi }{\hslash {\mathcal{A}}^2}\sum _{{\mathbf{k}}_1,\mathbf{q}\ne 0}{\left|V_{\mathbf{k},{\mathbf{k}}_1;{\mathbf{k}}_1-\mathbf{q},\mathbf{k}+\mathbf{q}}^{\mathrm{hh}}\right|}^2\left[{\mathcal{L}}_0({\epsilon }_{-({\mathbf{k}}_1-\mathbf{q})}^{\mathrm{h}}+{\epsilon }_{-(\mathbf{k}+\mathbf{q})}^{\mathrm{h}}-{\epsilon }_{-{\mathbf{k}}_1}^{\mathrm{h}}-{\epsilon }_{-\mathbf{k}}^{\mathrm{h}},{\mathsf{\Gamma }}_{\mathrm{h}})\right.\hfill \\ \hfill & \times & \left.\left\{f_{-({\mathbf{k}}_1-\mathbf{q})}^{\mathrm{h}}{f}_{-(\mathbf{k}+\mathbf{q})}^{\mathrm{h}}(1-f_{-{\mathbf{k}}_1}^{\mathrm{h}})+(1-f_{-({\mathbf{k}}_1-\mathbf{q})}^{\mathrm{h}})(1-f_{-(\mathbf{k}+\mathbf{q})}^{\mathrm{h}})f_{-{\mathbf{k}}_1}^{\mathrm{h}}\right\}\right]\hfill \\ \hfill & +& \frac{8\pi }{\hslash {\mathcal{A}}^2}\sum _{{\mathbf{k}}_1,\mathbf{q}\ne 0}{\left|V_{\mathbf{k},{\mathbf{k}}_1;{\mathbf{k}}_1-\mathbf{q},\mathbf{k}-\mathbf{q}}^{\mathrm{eh}}\right|}^2\left[{\mathcal{L}}_0({\epsilon }_{\mathbf{k}-\mathbf{q}}^{\mathrm{e}}+{\epsilon }_{-({\mathbf{k}}_1-\mathbf{q})}^{\mathrm{h}}-{\epsilon }_{\mathbf{k}}^{\mathrm{e}}-{\epsilon }_{-{\mathbf{k}}_1}^{\mathrm{h}},{\mathsf{\Gamma }}_{\mathrm{eh}})\right.\hfill \\ \hfill & \times & \left.\left\{f_{-({\mathbf{k}}_1-\mathbf{q})}^{\mathrm{h}}{f}_{\mathbf{k}-\mathbf{q}}^{\mathrm{e}}(1-f_{-{\mathbf{k}}_1}^{\mathrm{h}})+(1-f_{-({\mathbf{k}}_1-\mathbf{q})}^{\mathrm{h}})(1-f_{\mathbf{k}-\mathbf{q}}^{\mathrm{e}})f_{-{\mathbf{k}}_1}^{\mathrm{h}}\right\}\right].\hfill \end{array}$$

(40)

Computationally, the $\pi$-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\times 2$ Hamiltonian matrix are given by [58,59]

$${\epsilon }_{\mathbf{k}}^{\mathrm{e},\mathrm{h}}=\pm \frac{3}{2}{\gamma }_{0}bk\equiv \hslash v_{\mathrm{F}}k,$$

(41)

where ${\gamma }_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

$${\varphi }_{p_z,\mathbf{k}}^{(\mathrm{c},\mathrm{v})}\left(\mathit{r}\right)=\frac{1}{\sqrt{2}}\left[U_{\mathbf{k}}^{\left(1\right)}\left(\mathit{r}\right)\mp {\mathtt{e}}^{i{\theta }_{\mathbf{k}}}{U}_{\mathbf{k}}^{\left(2\right)}\left(\mathit{r}\right)\right],$$

(42)

where as shown in Equation (20), $U_{\mathbf{k}}^{\left(1\right)}\left(\mathit{r}\right)$ and $U_{\mathbf{k}}^{\left(2\right)}\left(\mathit{r}\right)$ are two sublattice Bloch functions built from the superposition of the periodic $2p_z$ orbitals, Ref. [59] and ${\theta }_{\mathbf{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]

$${\psi }_{p_z}\left(\mathit{r}\right)=C_{0}rcos\theta {\mathtt{e}}^{-Z^{\ast }r/2a_0},$$

(43)

where $C_0$ is a normalization factor, $a_0$ the Bohr radius, and an effective nucleus charge number $Z^{\ast }$ is $3.18$.

In particular, the structure factor introduced in Equation (38) can be calculated explicitly as

$$\begin{array}{ccc}\hfill {\mathcal{F}}_{{\mathbf{k}}^{\prime },\mathbf{k}}^{\left(s\right)}\left(\mathit{q}\right)& =& {\delta }_{{\mathbf{k}}^{\prime },\mathbf{k}+\mathbf{q}}\int d^2\mathit{r}{\left[{\mathsf{\Phi }}_{\mathbf{k}+\mathbf{q}}^{\left(s\right)}\left(\mathit{r}\right)\right]}^{\ast }exp\left(i\mathit{q}·\mathit{r}\right){\mathsf{\Phi }}_{\mathbf{k}}^{\left(s\right)}\left(\mathit{r}\right)\hfill \\ & \approx & {\delta }_{{\mathbf{k}}^{\prime },\mathbf{k}+\mathbf{q}}\langle {\varphi }_{p_z,\mathbf{k}+\mathbf{q}}^{\left(s\right)}\left(\mathit{r}\right)|exp(i\mathit{q}·\mathit{r})|{\varphi }_{p_z,\mathbf{k}}^{\left(s\right)}\left(\mathit{r}\right)\rangle ,\hfill \end{array}$$

(44)

where $s=c,v$ for Bloch wave function. Moreover, the Bloch-function structure factor in Equation (44) takes the form [60]

$$\begin{array}{ccc}\hfill & & \langle {\varphi }_{p_z,\mathbf{k}+\mathbf{q}}^{\left(s\right)}\left(\mathit{r}\right)|exp(i\mathit{q}·\mathit{r})|{\varphi }_{p_z,\mathbf{k}}^{\left(s\right)}\left(\mathit{r}\right)\rangle \hfill \\ & =& \frac{1}{N_A+N_B}\sum _{\mathbf{R}={\mathbf{R}}_A,{\mathbf{R}}_B}\langle {\psi }_{p_z}(\mathit{r}-\mathit{R})|exp\left[i\mathit{q}·(\mathit{r}-\mathit{R})\right]|{\psi }_{p_z}(\mathit{r}-\mathit{R})\rangle \frac{1}{2}\left[1\pm \frac{\gamma (\mathit{k}+\mathit{q}){\gamma }^{\ast }\left(\mathit{k}\right)}{\left|\gamma (\mathit{k}+\mathit{q})\gamma \left(\mathit{k}\right)\right|}\right],\hfill \end{array}$$

(45)

where tight-binding function ${\psi }_{p_z}\left(\mathit{r}\right)$ is given by Equation (43), and the signs (±) correspond to conduction (+) and valence (−) bands, respectively [59].

For intrinsic graphene, we have chemical potential ${\mu }_{\mathrm{e}}={\mu }_{\mathrm{h}}=0$ [61]. However, there is still a finite intrinsic areal density $n_{\mathrm{i}}\approx (\pi /6){(k_{\mathrm{B}}T/\hslash v_{\mathrm{F}})}^2$ due to thermal excitation of electrons and holes at finite temperatures T. In fact, we find $f_{\mathbf{k}}^{\mathrm{e}}=f_{\mathbf{k}}^{\mathrm{h}}=1/2$ at the K valley or $k=0$. Here, the calculated CDD rates from Equations (34) and (40), respectively, for electrons ${\mathsf{\Delta }}_{\mathrm{e}}\left(\mathit{k}\right)$ and holes ${\mathsf{\Delta }}_{\mathrm{h}}\left(\mathit{k}\right)$ are presented in Figure 9a at $T=77$ K and in Figure 9b at $T=300$K. Since $f_{\mathbf{k}}^{\mathrm{e},\mathrm{h}}\sim exp(-{\epsilon }_{\mathbf{k}}^{\mathrm{e},\mathrm{h}}/k_{\mathrm{B}}T)$ as ${\epsilon }_{\mathbf{k}}^{\mathrm{e},\mathrm{h}}\gg k_{\mathrm{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 ${\epsilon }_{\mathbf{k}}^{\mathrm{e},\mathrm{h}}$ around $k=0$, as seen in Figure 6.

The Coulomb diagonal-dephasing rates ${\mathsf{\Delta }}_{\mathrm{e},\mathrm{h}}\left(\mathit{k}\right)$ 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 ${\mathsf{\Lambda }}_{\mathrm{e},\mathrm{h}}(\mathit{k},\mathit{q})$ reveals deformations of induced optical-polarization waves with different wave vectors $\mathit{k}+\mathit{q}$ [8].

Considering the fact that major occupations of electrons and holes are accumulated around $k=0$, we have $f_{\mathbf{k}}^{\mathrm{e},\mathrm{h}}\approx 0$ only if k is large. As a result, we find from Equation (34) that $f_{{\mathbf{k}}_1-\mathbf{q}}^{\mathrm{e}}{f}_{\mathbf{k}+\mathbf{q}}^{\mathrm{e}}(1-f_{{\mathbf{k}}_1}^{\mathrm{e}})\ll 1$ at $k=0$ since we require $1-f_{{\mathbf{k}}_1}^{\mathrm{e}}\approx 1$ for large $k_1$, $f_{\mathbf{k}+\mathbf{q}}^{\mathrm{e}}\approx 1$ for small q, and $f_{{\mathbf{k}}_1-\mathbf{q}}^{\mathrm{e}}\approx 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_{{\mathbf{k}}_1-\mathbf{q}}^{\mathrm{e}}{f}_{-(\mathbf{k}-\mathbf{q})}^{\mathrm{h}}(1-f_{{\mathbf{k}}_1}^{\mathrm{e}})\ll 1$. Combining these two facts together, we expect that a dip will occur at $k=0$ for the Coulomb diagonal-dephasing rate ${\mathsf{\Delta }}_{\mathrm{e}}\left(\mathit{k}\right)$, as seen in Figure 9a. Moreover, the observed anisotropic energy dispersion in Figure 6a along the K-M and K-$\mathsf{\Gamma }$ directions directly leads to a staircase-like feature in Figure 9a for both ${\mathsf{\Delta }}_{\mathrm{e}}\left(\mathit{k}\right)$ and ${\mathsf{\Delta }}_{\mathrm{h}}\left(\mathit{k}\right)$. 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 ($\propto 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 ${\mathsf{\Delta }}_{\mathrm{e}}\left(\mathit{k}\right)$ and ${\mathsf{\Delta }}_{\mathrm{h}}\left(\mathit{k}\right)$ 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 ${\mathsf{\Delta }}_{\mathrm{h}}\left(\mathit{k}\right)$ in comparison with that of ${\mathsf{\Delta }}_{\mathrm{e}}\left(\mathit{k}\right)$, as well as different dispersion features around the K valley for ${\mathsf{\Delta }}_{\mathrm{e}}\left(\mathit{k}\right)$ and ${\mathsf{\Delta }}_{\mathrm{h}}\left(\mathit{k}\right)$. 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.

6. 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 intraband-scattering of holes can be obtained in a similar way.

By using the detailed-balance condition, the microscopic energy-relaxation time ${\tau }_{\mathrm{rel}}\left(\mathit{k}\right)$ of electrons in the presence of randomly-distributed ionized impurities can be calculated according to [44]

$$\frac{1}{{\tau }_{\mathrm{rel}}\left(\mathit{k}\right)}={\mathcal{W}}_{\mathrm{in}}\left(\mathit{k}\right)+{\mathcal{W}}_{\mathrm{out}}\left(\mathit{k}\right),$$

(46)

where the scattering-in rate for electrons in the final $|\mathit{k}\rangle$ state is

$$\begin{array}{ccc}\hfill {\mathcal{W}}_{\mathrm{in}}\left(\mathit{k}\right)& =& \frac{4\pi n_{\mathrm{im}}}{\hslash \mathcal{A}}\sum _{\mathbf{q}\ne 0}\left\{{\left|U_{\mathbf{k},\mathbf{k}-\mathbf{q}}^{\mathrm{im}}\left(\mathit{q}\right)\right|}^2{f}_{\mathbf{k}-\mathbf{q}}^{\mathrm{e}}{\mathcal{L}}_0({\epsilon }_{\mathbf{k}}^{\mathrm{e}}-{\epsilon }_{\mathbf{k}-\mathbf{q}}^{\mathrm{e}},{\mathsf{\Gamma }}_{\mathrm{e}})\right.\hfill \\ & +& \left.{\left|U_{\mathbf{k},\mathbf{k}+\mathbf{q}}^{\mathrm{im}}\left(\mathit{q}\right)\right|}^2{f}_{\mathbf{k}+\mathbf{q}}^{\mathrm{e}}{\mathcal{L}}_0({\epsilon }_{\mathbf{k}}^{\mathrm{e}}-{\epsilon }_{\mathbf{k}+\mathbf{q}}^{\mathrm{e}},{\mathsf{\Gamma }}_{\mathrm{e}})\right\},\hfill \end{array}$$

(47)

whereas the scattering-out rate for electrons in the initial $|\mathit{k}\rangle$ state takes the form

$$\begin{array}{ccc}\hfill {\mathcal{W}}_{\mathrm{out}}\left(\mathit{k}\right)& =& \frac{4\pi n_{\mathrm{im}}}{\hslash \mathcal{A}}\sum _{\mathbf{q}\ne 0}\left\{{\left|U_{\mathbf{k}+\mathbf{q},\mathbf{k}}^{\mathrm{im}}\left(\mathit{q}\right)\right|}^2(1-f_{\mathbf{k}+\mathbf{q}}^{\mathrm{e}}){\mathcal{L}}_0({\epsilon }_{\mathbf{k}+\mathbf{q}}^{\mathrm{e}}-{\epsilon }_{\mathbf{k}}^{\mathrm{e}},{\mathsf{\Gamma }}_{\mathrm{e}})\right.\hfill \\ & +& \left.{\left|U_{\mathbf{k}-\mathbf{q},\mathbf{k}}^{\mathrm{im}}\left(\mathit{q}\right)\right|}^2(1-f_{\mathbf{k}-\mathbf{q}}^{\mathrm{e}}){\mathcal{L}}_0({\epsilon }_{\mathbf{k}-\mathbf{q}}^{\mathrm{e}}-{\epsilon }_{\mathbf{k}}^{\mathrm{e}},{\mathsf{\Gamma }}_{\mathrm{e}})\right\}.\hfill \end{array}$$

(48)

Here, $n_{\mathrm{im}}$ represents the areal density of ionized impurity atoms in the crystal, and $|U_{\mathbf{k},{\mathbf{k}}^{\prime }}^{\mathrm{im}}{\left(\mathit{q}\right)|}^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

$${\left|U_{\mathbf{k},{\mathbf{k}}^{\prime }}^{\mathrm{im}}\left(\mathit{q}\right)\right|}^2=Z^{\ast 2}{\left|u_{\mathrm{c}}\left(q\right)\right|}^2{\left|\int d^2\mathit{r}{\left[{\mathsf{\Phi }}_{\mathbf{k}}^{\left(\mathrm{c}\right)}\left(\mathit{r}\right)\right]}^{\ast }exp(i\mathit{q}·\mathit{r}){\mathsf{\Phi }}_{{\mathbf{k}}^{\prime }}^{\left(\mathrm{c}\right)}\left(\mathit{r}\right)\right|}^2=Z^{\ast 2}{\left|u_{\mathrm{c}}\left(q\right)\right|}^2{\left|{\mathcal{F}}_{\mathbf{k},{\mathbf{k}}^{\prime }}^{\left(\mathrm{c}\right)}\left(\mathit{q}\right)\right|}^2,$$

(49)

where $Z^{\ast }$ is the charge number of ionized impurity atoms.

Substituting Equation (49) back into Equation (47), we obtain

$$\begin{array}{ccc}\hfill {\mathcal{W}}_{\mathbf{k}}^{\mathrm{in}}& =& \frac{4\pi n_{\mathrm{im}}{Z}^{\ast 2}}{\hslash \mathcal{A}}\sum _{\mathbf{q}\ne 0}\left\{f_{\mathbf{k}-\mathbf{q}}^{\mathrm{e}}{\mathcal{L}}_0({\epsilon }_{\mathbf{k}}^{\mathrm{e}}-{\epsilon }_{\mathbf{k}-\mathbf{q}}^{\mathrm{e}},{\mathsf{\Gamma }}_{\mathrm{e}}){\left|u_{\mathrm{c}}\left(q\right)\right|}^2{\left|{\mathcal{F}}_{\mathbf{k},\mathbf{k}-\mathbf{q}}^{\left(\mathrm{c}\right)}\left(\mathit{q}\right)\right|}^2\right.\hfill \\ & +& \left.f_{\mathbf{k}+\mathbf{q}}^{\mathrm{e}}{\mathcal{L}}_0({\epsilon }_{\mathbf{k}}^{\mathrm{e}}-{\epsilon }_{\mathbf{k}+\mathbf{q}}^{\mathrm{e}},{\mathsf{\Gamma }}_{\mathrm{e}}){\left|u_{\mathrm{c}}\left(q\right)\right|}^2{\left|{\mathcal{F}}_{\mathbf{k},\mathbf{k}+\mathbf{q}}^{\left(\mathrm{c}\right)}\left(\mathit{q}\right)\right|}^2\right\},\hfill \end{array}$$

(50)

$$\begin{array}{ccc}\hfill {\mathcal{W}}_{\mathbf{k}}^{\mathrm{out}}& =& \frac{4\pi n_{\mathrm{im}}{Z}^{\ast 2}}{\hslash \mathcal{A}}\sum _{\mathbf{q}\ne 0}\left\{(1-f_{\mathbf{k}+\mathbf{q}}^{\mathrm{e}}){\mathcal{L}}_0({\epsilon }_{\mathbf{k}+\mathbf{q}}^{\mathrm{e}}-{\epsilon }_{\mathbf{k}}^{\mathrm{e}},{\mathsf{\Gamma }}_{\mathrm{e}}){\left|u_{\mathrm{c}}\left(q\right)\right|}^2{\left|{\mathcal{F}}_{\mathbf{k}+\mathbf{q},\mathbf{k}}^{\left(\mathrm{c}\right)}\left(\mathit{q}\right)\right|}^2\right.\hfill \\ & +& \left.(1-f_{\mathbf{k}-\mathbf{q}}^{\mathrm{e}}){\mathcal{L}}_0({\epsilon }_{\mathbf{k}-\mathbf{q}}^{\mathrm{e}}-{\epsilon }_{\mathbf{k}}^{\mathrm{e}},{\mathsf{\Gamma }}_{\mathrm{e}}){\left|u_{\mathrm{c}}\left(q\right)\right|}^2{\left|{\mathcal{F}}_{\mathbf{k}-\mathbf{q},\mathbf{k}}^{\left(\mathrm{c}\right)}\left(\mathit{q}\right)\right|}^2\right\}.\hfill \end{array}$$

(51)

Using the inverse microscopic energy-relaxation time in Equation (46), we can further calculate the macroscopic thermally-averaged energy-relaxation time ${\tau }_{\mathrm{rel}}\left(T\right)$ as a function of temperature T, yielding [44]

$$\frac{1}{{\tau }_{\mathrm{rel}}\left(T\right)}=\frac{4}{n_{\mathrm{e}}\mathcal{A}}\sum _{\mathbf{k}}\left[\frac{1}{{\tau }_{\mathrm{rel}}\left(\mathit{k}\right)}\right]f_{\mathbf{k}}^{\mathrm{e}}.$$

(52)

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_{\mathrm{F}}=\hslash v_{\mathrm{F}}\sqrt{\pi n_0}$ at low temperatures, Ref. [61] where $n_0$ represents the areal electron density from doping, i.e., $n_0=n_{\mathrm{im}}$ for completely ionized doping atoms. For low temperatures with $k_{\mathrm{B}}T\ll E_{\mathrm{F}}$, we have $f_{\mathbf{k}}^{\mathrm{e}}=\mathsf{\Theta }(E_{\mathrm{F}}-{\epsilon }_{\mathbf{k}}^{\mathrm{e}})$ or $f_{\mathbf{k}}^{\mathrm{e}}=\mathsf{\Theta }(k_{\mathrm{F}}-k)$, where $\mathsf{\Theta }\left(x\right)$ is a unity step function and $k_{\mathrm{F}}=\sqrt{\pi n_0}$ is the Fermi wave number.

Physically, the Coulomb diagonal-dephasing rates $\mathsf{\Gamma }\left(\mathit{k}\right)={\mathsf{\Delta }}_{\mathrm{e}}\left(\mathit{k}\right)+{\mathsf{\Delta }}_{\mathrm{h}}\left(\mathit{k}\right)$ 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/{\tau }_{\mathrm{rel}}\left(T\right)$, 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 Figure 9 and Figure 10, respectively, represent two fundamentally different microscopic physics mechanisms.

As seen from Figure 10, we find the electron energy-relaxation rate $1/{\tau }_{\mathrm{rel}}\left(T\right)$ reduces with increasing temperature T due to enhanced screening effect on Coulomb interaction $u_{\mathrm{c}}\left(q\right)$ between two electrons or the rising of $q_0$ in Equation (37) with T, which implies that we have to wait a longer time ${\tau }_{\mathrm{rel}}\left(T\right)$ for our system returning to its initial thermal-equilibrium state at an elevated temperature. Furthermore, using the second-order Boltzmann moment equation [44], we would emphasize that this average energy-relaxation time ${\tau }_{\mathrm{rel}}\left(T\right)$, 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.

7. 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 first-principles 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.