Band dispersion of graphene with structural defects

We study the band dispersion of graphene with randomly distributed structural defects using two complementary methods, exact diagonalization of the tight-binding Hamiltonian and implementing a self-consistent T matrix approximation. We identify three distinct types of impurities resulting in qualitatively different spectra in the vicinity of the Dirac point. First, resonant impurities, such as vacancies or 585 defects, lead to stretching of the spectrum at the Dirac point with a finite density of localized states. This type of spectrum has been observed in epitaxial graphene by photoemission spectroscopy and discussed extensively in the literature. Second, nonresonant (weak) impurities, such as paired vacancies or Stone-Wales defects, do not stretch the spectrum but provide a line broadening that increases with energy. Finally, disorder that breaks sublattice symmetry, such as vacancies placed in only one sublattice, open a gap around the Dirac point and create an impurity band in the middle of this gap. We find good agreement between the results of the two methods and also with the experimentally measured spectra.

Graphene presents high potential for providing the next generation of electronic materials due to its strictly twodimensional character as well as its high electron mobility. It has demonstrated high design exibility, such as doping by atoms or molecules, e cient decoupling from an underlying substrate, or high tensile strength for exible electronics [1][2][3]. Several theoretical proposals as well as experiments are concerned with enhancing the spin-orbit coupling to open a band gap or inducing spin-spli ing [4,5]. Even a possible transition to a superconducting state has been proposed [6]. However, one of the most important prerequisites for graphene to become a base material for future electronics concerns the opening of a band gap, which has not been successfully demonstrated so far. e linear crossing of the bands near the Dirac point is protected by symmetry, because the two subla ices are equivalent. In order to open a band gap, this subla ice symmetry has to be broken.
Angle-resolved photoemission spectroscopy (ARPES) is the most direct method to probe the electronic structure experimentally. Numerous studies have examined the band structure of graphene near the Dirac point [7][8][9][10][11][12][13][14][15][16]. Several of these studies observe an elongated region near Dirac point as if the two touching cones are pulled apart, stretched but without tearing apart. Such occurrences have been discussed extensively in literature and depending on the speci c environment of the graphene, were a ributed to imperfections in the graphene, interactions with the substrate, or the opening of a band gap [16][17][18][19][20][21]. One observation common to all instances of the stretched Dirac point is the residual spectral weight that is still present at lowest energies. It needs to be understood in more detail in order to judge, if and under what conditions it can be referred to as an actual gap.
In this le er, we present a real-space tight-binding calculation modeling di erent kinds of structural defects randomly distributed over a graphene sample. We show that in all instances, except the case of vacancies placed in a single subla ice, there is no band gap opening near the Dirac point. e spectrum near the Dirac point is either almost unchanged or exhibits a stretched Dirac point with broadened states, which resemble experimentally observed band dispersion. Complementing our tight-binding model with a self-consistent T matrix approximation (SCTMA) calculation [22,23], we show that point defects in graphene are either a resonant or nonresonant type [24,25]. It is resonant defects that produce a dispersion with broadened states near the Dirac point resembling the results of the tight-binding calculation as well as experimental ndings. Due to the remarkable similarity between the results of the SCTMA, numerical simulations of the tight-binding model, and the experiment, we conclude that the broadened states measured in epitaxial graphene must at least in part be caused by speci c resonant defects in graphene. On the other hand, nonresonant defects are shown to weakly a ect the dispersion near the Dirac point, which is also con rmed by the tight-binding model calculation. ese ndings lead to the conclusion that it is close to impossible to open a band gap in graphene by virtue of defects only. We employed a second nearest-neighbor real-space tightbinding model with the following Hamiltonian: Here the indices i and j label individual atoms with one p zorbital per atom. e parameters t and t are hopping amplitudes between rst and second nearest neighbors, respectively.
e sums with single and double angular brackets run over rst and second nearest neighbors, respectively. We have typically built the Hamiltonian for a rectangular supercell of 160,000 carbon atoms with periodic boundary conditions and add randomly distributed structural defects with a given concentration. e Hamiltonian is then exactly diagonalized numerically and the resulting real-space wave functions are converted to momentum-space. is yields the complete set of eigenenergies E n and corresponding eigenfunctions ψ n (p). e spectral weight function is then computed arXiv:1811.00087v1 [cond-mat.mtrl-sci] 31 Oct 2018 e code was programmed in Matlab using the built-in routines for calculating eigenvalues and eigenvectors as well as Fourier transforms. e defects were introduced by suppressing the hopping between particular randomly chosen la ice sites. e atom positions near the defects have not been relaxed, except for the Stone-Wales (SW) defect, which involves repositioning of two carbon atoms [26]. e concentration of defects, n imp , is the ratio of the number of carbon atoms taken out or displaced to the total number of carbon atoms in the la ice. e large number of la ice sites used in the calculation is necessary for statistical reasons to increase the number of eigenstates in the region of low density of states (DOS) near the Dirac point as well as to have a large enough number of randomly distributed defects in the supercell. For the calculations presented in this le er the values of the parameters are t = 3.033 eV and t = 0.2 eV [27].
In the alternative SCTMA calculation, we start with the exact Green's function on a honeycomb la ice at zero energy. Only nearest-neighbor hopping terms are taken into account with t = 3.033 eV. Any individual structural defect considered in our calculation perturbs at most six neighboring sites of the la ice, hence we describe it with an exact 6 × 6 T matrix. We then convert this exact zero-energy T matrix to the basis of 4-component spinors governed by the continuous low-energy massless Dirac Hamiltonian with two valleys [28]. e converted T matrix is averaged over all possible positions and orientations of the defect and a non-zero energy is introduced as a perturbation. e averaged Green's functions of graphene with a nite concentration of defects, n imp , acquire the self-energy Σ = (n imp /A) T (E − Σ) . Here we have also included the same self energy in the argument of the impurity T matrix thus introducing the self-consistency equation. e spectral weight is calculated from the self energy as and is compared to the results of the tight-binding model and to the experimentally measured dispersion. Further details of the SCTMA calculation can be found in the Supplemental Material [23]. We consider ve di erent types of structural defects illustrated in Fig. 1. ey show the e ect of subla ice symmetry breaking and the qualitative di erence between resonant and nonresonant defects. e rst three types of impurities are vacancies that are either distributed in a single (A) subla ice, equally in both subla ices (A and B), or paired (whole AB unit cells removed). ese defects are shown in Fig. 1(a-c). While vacancies are not feasible in graphene, they provide a good model for adatoms a ached to individual la ice sites inducing a strong on-site potential [24,29]. Vacancies allow us to probe subla ice symmetry breaking and to demonstrate its e ect on the band dispersion [25,[28][29][30][31]. We also consider two other structural point defects, that have been experimentally observed in epitaxial graphene [32] and are shown in Fig. 1(d, e). e 585 defect, Fig. 1(d), is similar to the AB paired vacancy but with two reconstructed bonds. e SW defect is shown in Fig. 1(e) and involves a 90 • rotation of a bond between two adjacent atoms, along with the rearranging of the hopping terms around them.
Numerically computed spectra for distributions of vacancies are shown in Fig. 1(f-h). We see that vacancies placed in only one subla ice, Fig. 1(f), open a gap in the Dirac spectrum with an additional midgap band [31]. We will discuss the origin of this extra band below. Removing atoms randomly from both subla ices, as shown in Fig. 1(g), induces additional spectral weight in the vicinity of the Dirac point with momentum broadening that gets stronger closer to zero energy. At higher energies we see a band structure resembling "elongated" Dirac points. Finally, removing adjacent atoms, Fig. 1(h), does not noticeably a ect the spectrum apart from an energy dependent broadening of the line width. Among these three qualitatively di erent electronic spectra, only the rst one (vacancies placed in one subla ice, Fig. 1(a)) provides a true band gap, Fig. 1(f). e structural defects shown in Fig. 1(d, e) also demonstrate qualitatively di erent spectra similar to Fig. 1(g, h), respectively. We analyze them analytically in the framework of the SCTMA; see Supplemental Material [23]. For the 585 defect we nd the following equation for the self energy: is form is the result of a divergent zero-energy T matrix [22,23]. Such defects are known as resonant. For the case of the 585 defect β = 12.5 eV 2 . In Fig. 2 we compare the spectral weight obtained from direct numerical diagonalization of the disordered la ice model (panels (a-c)) with the solution of the self-consistency Eq (4) (panels (d-f)). For both calculations, the resulting structure near the Dirac point is very similar to experimentally measured graphene dispersion. Stretching of the spectrum near the Dirac point can be estimated from Eq. (4) as where c is a ing parameter. We plot the size of the smeared region around the Dirac point in Fig. 2(g) along with the t (c = 2.0712 ± 0.38).
Another type of structural defects (SW, Fig. 1(e)) belongs to the nonresonant case. In this case, the SCTMA provides the following self energy equation [23]: e nonresonant case is typically what is found for most point defects and for the SW defect α = 6.85. e band dispersion from the direct numerical diagonalization and from the SCTMA are shown in Fig. 3. ere is no apparent "elongated" Dirac point even at relatively high concentrations of impurities. Instead an energy-dependent line broadening gets stronger away from the Dirac point. is is in contrast to the relatively uniform broadening found in the band dispersion with resonant impurities. We see that both for resonant 585 defects and non-resonant SW defects, the results of direct diagonalization and SCTMA calculation show a remarkable agreement. e two forms of the self-energy, Eqs. (4) and (6), are the only two possibilities for point defects (as long as sublattice symmetry is preserved). In the resonant case, the zeroenergy T matrix diverges hence it can be approximated as T ∼ 1/E (up to a logarithmic factor).
is leads to the selfconsistency equation of the form in Eq (4). In the nonresonant case, small-energy expansion of the T matrix starts with a nonessential constant (it can be absorbed in the chemical potential) and a linear term leading to Eq. (6). Both resonant and nonresonant cases are in a good agreement with our direct numerical simulations of the tight-binding model. When the subla ice symmetry is broken, self energy is not a number anymore but rather an operator in the subla ice space. en more possibilities beyond Eqs. (4) and (6) emerge [25]. Vacancies distributed in only one subla ice, Fig. 1(a, f), is one possible illustration of this e ect.
We also consider a mixture of 585 and SW defects in a single sample to gain insight on how these defects interact and to have an approximation of realistically disordered graphene. In Fig. 4, we compare the band structure for an equal amount of 585 and SW defects obtained from direct diagonalization (panel (a)) and the SCTMA calculation is behavior is typical for nonresonant impurities.
(panel (b)) with the ARPES data of epitaxial graphene (panel (c)) [16][17][18][19][20][21]. e average DOS is shown in Fig. 4(d, e) for the tight-binding model and SCTMA, respectively. In the tightbinding model we implement an equal number of 585 and SW defects in the la ice while for the SCTMA calculation the self-consistency equation contains the sum of the right-hand sides of Eqs. (4) and (6). e DOS is a momentum integral of the previously calculated spectral weight. e results in Fig. 4(a, b) not only show a striking resemblance to each other but also match the "elongation" found in experimentally measured epitaxial graphene dispersions as seen in Fig. 4(c). e DOS shows a nearly constant region near the Dirac point while further away it matches the linear dispersion of states expected for graphene.
is nite DOS near the Dirac point is dominated by the e ect from resonant defects whose T matrix diverges at zero energy.
e DOSs also show an excellent agreement between the SCTMA and the direct tight-binding model.
From our calculations, we can make several conclusions and comment on some new insight. Firstly, the only possibility for a gap in the spectrum is when the subla ice symme- try is broken, for example in the case of vacancies in a single subla ice, Fig. 1(f). Mathematically, this can be understood by representing the tight-binding Hamiltonian in the matrix form as [31] Here, the block h contains hopping terms from one sublattice to the other and h † de nes the opposite hopping. Second nearest neighbor hopping is temporarily disregarded. Removing atoms from only one subla ice causes h and h † to be non-square blocks. e matrix (7) has a number of zero eigenvalues equal (or larger [33]) to the di erence of the number of A and B sites. Hence the DOS exhibits a delta-peak at zero energy with a gap opening around it. is gap is a result of statistical repulsion between zero and non-zero eigenstates of the random Hamiltonian matrix [30,31]. When the second nearest neighbor hopping is taken into account, the zeroenergy eigenstates rearrange into a dispersive midgap band found in Fig. 1(f). ese results suggest that selective removal of (or chemical bonding to) carbon atoms from a given subla ice may be the only way to realize gapped graphene experimentally. e second insight is about the origin of "elongated" Dirac point that has been measured experimentally in epitaxial graphene [16][17][18][19][20][21]. e SCTMA shows that a point defect in graphene can either be resonant or nonresonant, where the geometry of a defect a ects only the parameters α or β. Due to the remarkable similarity between experimentally measured "elongated" Dirac point (Fig. 4(c)), the tight-binding result ( Fig. 4(a)) and the SCTMA result ( Fig. 4(b)), we conclude that this "elongation" is caused, at least in part, by resonant defects. We have shown that 585 defects are resonant and provide an apparent stretching of the Dirac point. At the same time, it is known that 585 defects are common in epitaxial graphene [32]. We thus conclude that the "elongation" of the Dirac point observed in many graphene samples is the result of resonant defects and 585 defects in particular. Regarding the prospective band gap in epitaxial graphene, our calculations explicitly show that there are states in the Dirac point region, leading us to conclude that the "elongation" of the Dirac point can not be considered as a gap. Furthermore, any concentration of resonant defects will increase the number of states at low energies. Nevertheless, the stretching of the spectrum near the Dirac point creates an energy range were electrons are localized in real space. is phenomenon may be used to open a mobility gap around the Dirac point in graphene and realize an insulating state [34].
In summary, we implemented a simple real-space tightbinding model that allowed us to calculate the band dispersion of graphene with defects. We found that a band gap can only be induced when the subla ice symmetry is broken as shown in Fig. 1(f). Looking at more realistic defects, we found that the 585 defect creates an "elongated" Dirac point (Fig. 2) similar to those found in experimentally measured spectra of epitaxial graphene. Graphene with SW defects has a qualitatively di erent spectrum without apparent stretching (Fig. 3). We conclude that the experimentally observed "elongated" Dirac point is due to the 585 defects and the length of stretching scales as the square root of their concentration. We have also developed a SCTMA theory [23] which allowed us to classify point defects as either resonant or nonresonant and showed a remarkable agreement with the direct la ice calculations.
e SCTMA model provides further insights into the nature of the broadened states measured in epitaxial graphene, showing that the "elongated" region can not be considered a spectral gap. At the same time, disorder can lead to localization of the states near the Dirac point and hence to a mobility gap [34].
Rev. B 94, 064204 (2016). In this Supplemental Material we provide a detailed derivation of the self-consistent T matrix approximation for structural defects in graphene. We derive exact analytical zero-energy T matrices for Stone-Wales and 585 defects within the tight-binding model. ese results are then translated to the e ective low-energy description of graphene with the Dirac Hamiltonian. Finally, averaging with respect to positions and orientations of the defects is performed within the self-consistent approximation.

A. Tight-binding model
We describe electrons in graphene within the nearest neighbor tight-binding model (second nearest neighbor hopping, cf. Eq. (1) of the main text, will be neglected throughout this calculation). Wave functions Ψ(r) are de ned on the honeycomb la ice with the la ice spacing a = 2.46Å; the spatial argument r takes the corresponding discrete values. We assume the hopping amplitude of t = 3.033 eV as in the main text of the paper. e tight-binding Hamiltonian acts according tô Here r refers to la ice sites, ζ r = ±1 is a sign function distinguishing A and B subla ices, and vectors δ s refer to the three nearest neighbors e angle α de nes orientation of the crystal with respect to coordinate axes. e two vectors of the elementary la ice translations are δ 1 − δ 0 and δ 2 − δ 0 .
Zero-energy Green's function is simply the inverse of the Hamiltonian operator. It has nonzero matrix elements only between sites from di erent subla ices. We construct the Green's function in momentum representation and then transform it to real-space as e argument of this function is the vector connecting a site in the A subla ice with a site in the B subla ice displaced by n 1,2 elementary translations in the directions δ 1,2 − δ 0 . e rst few values of the Green function are displayed in Fig. S1.
We will consider the Stone-Wales defect and the 585 defect that perturb at most six neighboring sites of the la ice (shaded region in Fig. S1). We select a basis of these six sites to be e o set position r 0 de nes location of the impurity and refers to some site in the A subla ice. Between these six points, the la ice Green's function takes the valueŝ ,

B. Dirac Hamiltonian
E ective description of the low-energy electrons in graphene is provided by the massless Dirac Hamiltonian. In the valleysymmetric form it can be wri en as At short distances ϵr this function has the asymptotics G(iϵ, r) ≈ −iσ r 2π r 2 + iϵ 2π 2 2 log where γ ≈ 0.577 is the Euler-Mascheroni constant.
C. Relation between tight-binding and Dirac description e four-component wave function |Φ(r) governed by the Dirac Hamiltonian is related to the la ice wave function Ψ(r) in the following way (see Ref. S1) Here A is the area of the unit cell. (S10) e diagonal unitary matrix U encodes dependence on the impurity position r 0 and the la ice orientation α.

D. Impurity T matrix
A single impurity is described by the perturbationV to the tight-binding Hamiltonian. At zero energy, it corresponds to the T matrix de ned on the la ice ast (S11) In the Dirac language, the T matrix becomes Here we have also introduced the notation T 0 for the T matrix of an impurity placed at r 0 = 0 with the orientation α = 0. e general T matrix di ers from T 0 by the diagonal unitary rotation U only. At a nite Matsubara energy E = iϵ, we can express the T matrix using the following identities: Here ∆G(iϵ) denotes the di erence of two Green's functions of the Dirac Hamiltonian (S8) taken at coincident points.
∆G(iϵ) = lim r→0 G(iϵ, r) − G(0, r) ≈ iϵ log(ϵ/t) 2π 2 2 . (S14) Here we have used Eq. (S8) and replace r → a in the argument of the logarithm at short distances. E. Self energy for weak impurities: Stone-Wales defect We rst apply the above formalism to graphene with Stone-Wales defects. An isolated impurity is equivalently described by one of the two perturbation operatorŝ (S15) From the perturbation matrix we construct the la ice T matrixt at zero-energy using Eq. (S11) and translate it into the Dirac language applying Eq. (S12). For the moment we disregard the phase factors U † and U . e resulting T matrix has only two non-zero eigenvalues and can be represented as