Point Defects in Twisted Bilayer Graphene: A Density Functional Theory Study

Point Defects in Twisted Bilayer Graphene: A Density Functional Theory Study

Abstract

We have used ab initio density functional theory, incorporating van der Waals corrections, to study twisted bilayer graphene (TBLG) where Stone-Wales defects or monovacancies are introduced in one of the layers. We compare these results to those for defects in single layer graphene or Bernal stacked graphene. The energetics of defect formation is not very sensitive to the stacking of the layers or the specific site at which the defect is created, suggesting a weak interlayer coupling. However signatures of the interlayer coupling are manifested clearly in the electronic band structure. For the “” Stone Wales defect in TBLG, we observe two Dirac cones that are shifted in both momentum space and energy. This up/down shift in energy results from the combined effect of a charge transfer between the two graphene layers, and a chemical interaction between the layers, which mimics the effects of a transverse electric field. Charge density plots show that states near the Dirac points have significant admixture between the two layers. For Stone Wales defects at other sites in TBLG, this basic structure is modified by the creation of mini gaps at energy crossings. For a monovacancy, the Dirac cone of the pristine layer is shifted up in energy by eV due to a combination of the requirements of the equilibration of Fermi energy between the two layers with different numbers of electrons, charge transfer, and chemical interactions. Both kinds of defects increase the density of states at the Fermi level. The monovacancy also results in spin polarization, with magnetic moments on the defect of 1.2 – 1.8 .

pacs:
73.22.Pr, 31.15.A-, 61.72.-y, 71.15.Mb

I Introduction

The mechanical and electronic properties of graphene make it an attractive candidate for use in the electronics industry, and indeed it has even been speculated that it could one day replace silicon as the primary material used.(1); (2) In addition, it exhibits several exotic properties, such as a room temperature quantum Hall effect.(1) Many of the interesting electronic properties of single-layer graphene arise from the presence in the band structure of Dirac cones at the K points in the Brillouin zone that touch at the Fermi level , as a result of which electrons behave like massless Dirac fermions.(3)

When graphite or multilayer graphene is assembled by arranging single graphene layers in the conventional Bernal stacking, also referred to as AB-stacking, where successive layers are displaced by both a vertical and a lateral shift, many of the unique properties of single-layer graphene, such as the Dirac cones, disappear. However, in recent years there has been the appealing discovery that if instead a few layers of graphene are grown so that successive layers are aligned with a twist with respect to each other, there then appears to be an effective electronic decoupling between layers, so that the Dirac cones are maintained, at least for large angles of twist.(4); (5); (6) Twisted bilayer graphene has also been shown to have other intriguing properties, e.g., for small angles of twist, there is a reduction in the Fermi velocity, and a localization of electrons, and one can bring van Hove singularities in the electronic density of states very close to the Fermi level.(7) This in turn raises the possibility of seeing other interesting phenomena such as superconductivity. Other fascinating results on TBLG include the demonstration of Hofstadter’s Butterfly in the energy spectrum in a magnetic field,(8) and neutrino-like oscillations as a result of coupling of the Dirac cones of the two rotated layers.(9) Graphene layers with a twist can now be grown by a variety of methods, which allow one to access a range of twist angles. (10); (11); (12) However, the nature and extent of the interlayer coupling, and the consequences thereof, continue to be debated.(13)

The question then arises: in what way are these properties altered upon the introduction of defects in TBLG? There has always been a lot of interest in defects in graphene and other low-dimensional carbon materials.(14) Much of this interest stems from concerns about how the presence of such defects will affect the functioning of these materials in applications. Defects such as monovacancies, divacancies, interstitials and Stone-Wales defects(15) are known to affect the mechanical, electronic and magnetic properties of carbon materials. There have, for example, been several studies to check to what extent the presence of defects affects electrical conductance.(16); (17); (18) However, it has also been appreciated that it might be possible to purposely design defective structures with a view toward creating certain functionalities.(17); (19) Among the possible advantages that the presence of defects in graphene can confer, it has been shown that they can induce magnetic moments,(20); (21); (22); (23); (24); (25) and improve gas sensing properties.(26) It has also been suggested that the introduction of defects can be used for band-gap engineering.(27); (17); (28) Defects can be deliberately created by, e.g., irradiation with electrons(29); (30); (31); (32) or ions,(33); (34); (35); (36) oxidation,(37) hydrogenation,(38); (39); (40); (41) or fluorination.(42); (43)

In this theoretical study we perform density functional theory (DFT) calculations to study the properties of two kinds of point defects: Stone-Wales defects (SW) and monovacancies in twisted bilayer graphene. We are particularly interested in the question of how the electronic band structure in the neighborhood of ( i.e., the Dirac cones) are affected by the presence of defects.

Ii Computational Details

Our spin-polarized DFT calculations were performed using the PWscf package of the Quantum ESPRESSO distribution.(44) Ultrasoft pseudopotentials(45) were used to describe the interactions between the ion cores and valence electrons. A plane wave basis set was used, with kinetic energy cut-offs of Ry and Ry for the wavefunctions and charge densities, respectively. We have considered two different levels of approximation for the exchange-correlation interactions: (i) the local density approximation (LDA) in the Perdew-Zunger form, (46) and (ii) the generalized gradient approximation (GGA) in the Perdew-Burke-Ernzerhof (PBE) form.(47) Many of the previous DFT studies on TBLG made use of the LDA. However, given that the (weak) coupling between the two graphene layers is important in TBLG, and that there has been much debate about its nature and extent, we feel that it is desirable to have an accurate treatment of these interlayer interactions, where van der Waals interactions (which are absent in conventional DFT calculations) may be expected to play a crucial role. For this reason, along with the PBE exchange-correlation, we also incorporate the van der Waals interactions as a semi-empirical correction, using the “DFT-D2” treatment suggested by Grimme.(48); (49) Of the three kinds of theoretical treatments used in this paper (LDA, PBE and DFT-D2), we believe that those results using DFT-D2 should be regarded as being the most reliable, and most of our results have been obtained using this approach. However, we also present some results with the other two treatments for purposes of comparison; we believe that this is of interest since several authors still continue to use these other functionals, especially LDA, in order to treat such systems, and it is therefore worth examining their reliability.

System geometries are described in Section III.2 below. Along the direction (normal to the graphene sheets) periodic images were separated by a vacuum spacing of about 13 Å. The Brillouin zone was sampled with a Monkhorst-Pack mesh of k-points for the primitive unit cell of graphene, and proportionately equivalent meshes for larger supercells of graphene systems considered in this paper, when performing the self-consistent-field calculations. However, when performing the post-processing calculations, when an extremely accurate k-point sampling was needed in order to obtain a precise computation of the Fermi energy, much finer grids were used, and the grid spacing was further reduced in regions of the Brillouin zone (near the Dirac points) where it was found necessary to have a particularly dense sampling; the weights assigned to k-points were appropriately adjusted. Convergence was aided by making use of the Marzari-Vanderbilt smearing scheme, (50) with a small smearing width of 0.001 Ry. All the structures were relaxed using Hellmann-Feynman forces (51) until the force on each atom was less than 0.001 Ry/bohr along all the three directions. Simulated scanning tunneling microscopy (STM) images were obtained using the Tersoff-Hamann theory (52) for STM images, as incorporated in the “pp.x” post-processing tool supplied with Quantum-ESPRESSO.

Iii Twisted Bilayer Graphene

Modeling twisted bilayer graphene, i.e., a bilayer of graphene with a rotational stacking fault, requires the construction of a commensurate supercell with the necessary finiteness required for practical computation. A number of such commensurate supercells are possible depending on the angle of rotation between the two graphene layers. (53) The derivation of the commensuration condition and the construction of the commensurate unit cell for an arbitrary set of angles between the two graphene layers is described below. Though this commensuration condition and the resulting Brillouin zone have been described earlier in Ref. (53), we briefly summarize them here, in order to make clear the notation used in this paper, and to enable the reader to follow the analysis of results that is presented later in Section IV.2.2.

iii.1 Geometric Properties of Twisted Bilayer Graphene:

A TBLG system is created by placing two layers of graphene on top of each other, one of which is “unrotated” (U), and the other “rotated” (R). In order to describe the TBLG system, one must find a commensurate unit cell for both the layers.

Let and be the primitive lattice vectors of the (unrotated) graphene honeycomb lattice. A supercell for the unrotated layer can be constructed by the supercell lattice vectors: and , where and are integers. The primitive lattice vectors of the second graphene layer, which has been rotated by an angle , are given by: and , where is the rotation matrix. For such a rotated lattice, a supercell is given by the lattice vectors: and , where and are integers.

The condition for commensuration is given by:

 n1aR1+n2aR2=m1aU1+m2aU2 (1)

For a standard choice of graphene primitive lattice vectors and , where is the lattice constant, the above commensuration relation becomes:

 (n1n2)=(cosθ+1√3sinθ+2√3sinθ−2√3sinθcosθ−1√3sinθ)(m1m2) (2)

This maps one integer pair () to another (), with the constraint that . As pointed out in Ref. (53), this is a Diophantine problem, which has to be solved to find out integer pairs () and () that satisfy Eq. (2). For integer pairs satisfying the above commensuration condition, the twist angle is given by:

 θ=cos−1(2m1n1+m1n2+m2n1+2m2n22(m21+m22+m1m2)) (3)

The primitive lattice vectors of the twisted bilayer, and , are related to the unrotated primitive lattice vectors of the single layer, and , by:

 AB1=m1aU1+m2aU2, (4)

and

 AB2=−m2aU1+(m1+m2)aU2. (5)

The smallest possible such unit cell corresponds to , ; by Eqs. (2) and (3), this gives , , and . This is the angle of twist used for the twisted bilayer graphene considered in our study, its primitive unit cell contains 14 atoms in each layer. The unit cell vectors for this primitive unit cell are and .

In Fig. 1 we have shown the structure of TBLG with this twist angle, with the primitive unit cell of the twisted bilayer indicated by the black rhombus. Note that for this choice of twist angle, the upper (rotated) layer contains three symmetry-inequivalent types of carbon atoms, which we label , and (see Fig. 1).

However, most of our calculations have been performed using a larger unit cell , which has lattice vectors that are twice as long: and . This unit cell contains 56 atoms in each layer, i.e., 112 atoms in the bilayer unit cell; this corresponds to the choice , .

It is important to also have a knowledge of the structure of the reciprocal space of TBLG, in order to gain an understanding of its electronic properties. If and are the reciprocal lattice vectors corresponding to the unrotated layer, and and are the reciprocal lattice vectors corresponding to the rotated layer, the reciprocal lattice vectors of the TBLG are given by:

 GB1={m1+m2m21+m22+m1m2}gU1+{m2m21+m22+m1m2}gU2 (6)
 GB1=⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩37gU1+17gU2,for S1(6a)314gU1+114gU2,for S2(6b)

and

 GB2=−{m2m21+m22+m1m2}gU1+{m1m21+m22+m1m2}gU2 (7)
 GB2=⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩−17gU1+27gU2,for S1(7a)−114gU1+214gU2,for S2(7b)

Thus,

 gU1=m1GB1−m2GB2 (8)
 gU1=⎧⎪ ⎪⎨⎪ ⎪⎩2GB1−GB2,for S1(8a)4GB1−2GB2,for S2(8b)

and

 gU2=m2GB1+(m1+m2)GB2 (9)
 gU2=⎧⎪ ⎪⎨⎪ ⎪⎩GB1+3GB2,for S1(9a)2GB1+6GB2,for S2(9b)

The reciprocal lattice and the first Brillouin zone for the unrotated and rotated layers, and the twisted bilayer graphene, are shown in Fig. 2. The first Brillouin zone for the TBLG (shown by the black hexagon in Fig. 2) corresponds to a twist angle of , for the smallest commensurate supercell . A similar figure can be obtained when is used; this is provided in the Supplementary material for easy reference. Note that by the translational symmetry of the lattice, , and are identical; similarly, , and are identical. If, in addition, one has inversion symmetry or time-reversal symmetry, then all six K points are identical. However, these symmetries are valid at the corners of the Brillouin zone only, and need not necessarily hold in the interior of the Brillouin zone, i.e., one need not have three-fold or six-fold symmetry always present. In particular, the introduction of defects may lower the symmetry. The and points in the Brillouin zone of the unrotated layer can be folded onto points in the first Brillouin zone of the TBLG lattice:

 KU1=(m1−m2)KB1+m2(GB1+GB2) (10)
 KU1=⎧⎪ ⎪⎨⎪ ⎪⎩KB1+(GB1+GB2),for S1(10a)KB2+(2GB1+3GB2),for S2(10b)

and

 KU2=(m1−m2)KB2+m2GB1 (11)
 KU2=⎧⎪ ⎪⎨⎪ ⎪⎩KB2+GB1,for S1(11a)KB1+3GB1,for S2(11b)

Similarly for the rotated layer,

 KR1=(n1−n2)KB1+n2(GB1+GB2) (12)
 KR1=⎧⎪ ⎪⎨⎪ ⎪⎩KB1+GB2,for S1(12a)KB2+3GB2,for S2(12b)

and

 KR2=(n1−n2)KB2+n2GB1 (13)
 KR2=⎧⎪ ⎪⎨⎪ ⎪⎩KB2+(GB1+GB2),for S1(13a)KB1+(3GB1+2GB2),for S2(13b)

Thus, K folds onto , K or K, depending on whether () is , or , respectively. Similarly, for the rotated layer, K folds onto , K or K, depending on whether () is , or , respectively. For TBLG with the supercell , we therefore have the following mappings: K K and K K. Similarly, when performing calculations with the supercell , we have: K K and K K.

This tells us that when we use or to analyze the band-structure of the TBLG, the Dirac cones will continue to appear at the K points of the Brillouin zone of the supercell. Thus, these are the points whose vicinity we will focus on, to see whether the introduction of defects opens up a band gap and/or alters the linear dispersion relation.

iii.2 Systems Studied

We consider two types of point defects in TBLG, viz., the Stone-Wales defect and the monovacancy. We primarily work with the supercell which is times as large as the supercell . For both the Stone-Wales defect and the monovacancy, we perform calculations considering all possible inequivalent sites of one defect in the supercell. This corresponds to one defect in a supercell of 56 carbon atoms per layer, i.e., a defect density of 1.8%. For the monovacancy, we have in addition considered defect densities of 7.2% (corresponding to one vacancy in the supercell with 14 carbon atoms in each layer) and 0.8% (corresponding to one defect in supercell , which is times as large as , and contains 126 carbon atoms in each layer).

For comparison of various properties of the above-mentioned defects in TBLG, we also consider the same defects in a single layer of graphene (SLG), as well as in an AB-stacked bilayer graphene (we refer to this as AB-BLG). To facilitate comparison, i.e., maintain the same defect density and system symmetry, we use a similar supercell also for our calculations on SLG and AB-BLG.

Iv Results and Discussion

iv.1 Pristine Graphene

First, as benchmarks, we perform calculations on defect-free bilayer graphene, with both Bernal Stacking and twisted stacking. Our results for the energetics and geometry, for AB-BLG and TBLG, with both LDA and DFT-D2 approaches, as well as with the PBE alone, are shown in Table 1. While both LDA and DFT-D2 give values of the interlayer distance that are in reasonably good agreement with experiment, use of the PBE alone leads to a value that is far too large. However, the DFT-D2 leads to a result for the exfoliation energy , the energy required for exfoliating a graphene layer from the bilayer system, that is closer to experimental estimates than the LDA result; once again the PBE result is completely erroneous. All of these results are in agreement with previous theoretical work.(10); (58); (59) Also note that both LDA and DFT-D2 show that AB-BLG is lower in total energy than TBLG by about 2 meV/atom; however PBE alone leads to TBLG being more stable by the very small amount of 0.04 meV per atom. On going from AB-BLG to TBLG, there is a very slight increase in , which causes a small decrease in ; however, these changes are very small.

In Fig. 4 we plot the electronic band structure, as obtained using DFT-D2 for pristine (defect-free) twisted bilayer graphene, using the supercell, in the vicinity of the K point. In agreement with previous authors,(53) we obtain a conical dispersion, with the Dirac point lying exactly at K, and the Dirac crossing energy =. We have plotted the energy surfaces of the two topmost valence bands, we call these VB2 (blue in Fig. 4), VB1 (green), and the two lowest-lying conduction bands, CB1 (red) and CB2 (black). Note that VB1 and VB2 are almost, but not quite, degenerate, being shifted in energy by 9 meV; similarly CB1 and CB2 are displaced in energy by 9 meV. Note also that there is no band gap between VB1 and CB1. These results are similar to those obtained by previous authors who have performed calculations on TBLG.(53) Fig. 4 is useful in that it serves as a baseline to compare our other results to further below, when we will see how its features are modified upon the introduction of defects.

Figs. 2 (a)–(c) show simulated STM images of SLG, AB-BLG and TBLG; the input local densities of states for these calculations were obtained from DFT-D2 calculations. For a single layer of graphene, the STM image clearly shows the honeycomb like arrangement of carbon atoms. In the case of AB-stacked bilayer graphene, the image reflects the three-fold symmetry of Bernal stacking. These images are in excellent agreement with earlier reported experimental STM images. (60) In the case of TBLG, at first sight, the STM image might appear rather similar to that of the single layer of graphene, however closer examination of the STM image brings out the effect of the twisted lower layer. The relatively brighter spots in this image correspond to the positions where the atoms of both the layers lie directly atop each other, thus enhancing the STM signal.

iv.2 Stone-Wales Defect

A Stone-Wales (SW) defect is formed by rotating a single carbon-carbon bond in the graphene sheet by , resulting in a structure with a pair each of seven-membered and five-membered rings. This is known to be one of the most common defects in graphene-related systems and carbon nanotubes.(15) Since it has been reported from calculations on SLG (61) that an out-of-plane distortion of the carbon atoms in the vicinity of the SW defect further stabilizes this defect, we explicitly check if permitting such distortions stabilizes the SW defect in TBLG.

Structural properties

We create the SW defect in one of the layers of the TBLG, and compare the structural and electronic properties with those of SLG and AB-BLG with SW defects. Unlike SLG or AB-BLG, where only one type of site for a SW defect is possible, in TBLG, there can be many inequivalent defect-sites, with the number depending on the angle of twist between the two layers. For the case considered by us, where the two layers are rotated by a relative angle of 38.213, there can be four inequivalent geometries (see Fig. 5), corresponding to four different choices of the rotated bond; we name these -SW, -SW, -SW, and -SW, with the labeling convention denoting which pair of adjacent atoms is connected by the rotated bond (see Fig. 1).

For SLG, we find that a sinusoidal distortion pattern about the defect center results in the most stable geometry of the SW defect [see Figs. 6(a) and (b)], in agreement with earlier reports. (61) The amplitude of distortion, , (the difference between the largest upward displacement and the largest downward displacement) is as large as 1.2 Å [see Fig. 6(c)] for a SW defect in SLG within a supercell of size , which is in good agreement with earlier values of about 1 Å for supercells of similar size. In the case of the bilayers, we find that SW defects in both AB-BLG and TBLG are stabilized by a similar sinusoidal distortion pattern, however the amplitude of the distortion in the defective layer is reduced to about 0.7 Å; the presence of the second layer inhibits the distortion. In our study, we have also allowed the undefective layer to relax, it exhibits a smaller distortion amplitude of 0.2 Å. Fig. 6(c) shows the comparison of the net distortion of the graphene layers for the cases of single layer, AB-stacked bilayer, and twisted bilayer of graphene. Both the twisted and the AB-stacked bilayers show a very similar extent of distortion.

The energy of formation of the SW defect, as obtained with the various approaches, is given in Table 2. The values listed in parentheses in Table 2 give the formation energy for the SW defect when no out-of-plane distortion is permitted. Firstly for SLG, by comparing the results when this distortion is permitted with those when it is suppressed, we find that the sinusoidal distortion leads to a stabilization of about 220 meV/defect; this is in good agreement with earlier results.(61)

In the case of bilayers, the distortion results in a stabilization of the SW defect, as compared to the undistorted case, by an energy of about 40 meV in AB-BLG, and about 60 meV in TBLG. For TBLG, we find that the defect formation energy depends only very slightly on the position of the defect, this is because the interaction between the layers is small. Similarly, the energy to form a SW defect in TBLG is very similar to that in AB-stacked BLG. The small differences in defect formation energy (of the order of tens of meV) between the different kinds of Stone-Wales defects in TBLG are due to the differences in interlayer coupling when the SW defect is situated at different positions in the graphene sheet. Note that these differences disappear when the calculations are performed with the PBE alone, since in this case the interlayer coupling is described very poorly and becomes negligible. From the DFT-D2 calculation, we see that it costs an additional 140–180 meV to create a SW defect in TBLG than in SLG; the majority of this increase comes because the energy-lowering due to distortion along the -direction is hindered by the presence of the second layer. The same effect is reflected in our LDA results; once again, it is absent in the PBE results.

Band Structure

The band structure of AB-stacked bilayer graphene is characterized by parabolic bands that touch near the K points in the Brillouin zone, at the Fermi level . In contrast, for twisted bilayer graphene, there are Dirac cones, i.e., linear bands that touch at the K points, at (recall Fig. 4 above). Given the great interest in linear dispersion, as well as the potential importance of the band-gap engineering of graphene, we wish to know whether, upon the introduction of defects such as Stone-Wales defects and monovacancies, (i) the Dirac cones are maintained (ii) if not, in what way are they altered, i.e., are they shifted in energy (iii) whether gaps open up, and if so, at what energies. As we will show below, clear signatures of the importance of interlayer coupling emerge when looking at the band structure, as opposed to the energetics of defect formation.

We will examine the electronic band structure of SLG, AB-BLG and TBLG, containing a single SW defect in the supercell. No perceptible spin-polarization was found in any of the cases, and we therefore present non-spin-polarized results. Our structures break the reflection symmetry about the two axes passing through the SW defect. As a result, the irreducible Brillouin zone (IBZ) is half of the first Brillouin zone [shown by the shaded region in the TBLG Brillouin zone in Fig. 7]. Let us consider the high-symmetry points that lie in this IBZ. As mentioned in Section III above, by the translational symmetry of the lattice, K is identical to K; further, by making use of time-reversal symmetry, these are also identical to K. Thus, electronic eigenvalues should be identical at these three K points. On the other hand, no symmetries relate the three points M, M and M. However, as we will see below, it is not particularly useful, and perhaps even misleading, to plot the band structure for these systems by taking cuts along high-symmetry directions of the Brillouin zone, as is frequently done. Instead, we plot the surfaces of energy dispersion , in the vicinity of the point K, which is the region of especial interest. These results, for SW defects in SLG and AB-BLG, and the four distinct types of SW defects in TBLG, are plotted in Fig. 8.

Before we proceed to the band structure of TBLG, it is of interest to first consider the band structure of SLG with a SW defect, as it will help us interpret our results for TBLG with SW defects. The band-structure of SW defects in SLG has been the subject of a debate in the literature. While some authors reported that the presence of the SW defect opened up a gap,(28); (63) other authors have shown that this phenomenon depends on the size of the unit cell used;(64) in addition, there is a shift of the Dirac point away from the K point. The cause of these shifts in the Dirac point has been attributed to electron-phonon coupling;(65) note that the 90 bond rotation involved in the formation of SW defects can be viewed as a linear combination of -point phonon modes.(64) In Fig. a we see that the Dirac cone is preserved even after introducing the SW defect, but the Dirac point has been shifted in k-space [see also Fig. 7, where the position of the shifted Dirac point is indicated by the asterisk]. Note however that the Dirac crossing energy remains at the Fermi energy .

Next, we consider the band structure of a SW defect in AB-BLG [see Fig. b]. The higher valence band VB1, and the lower conduction band CB1, have been plotted; VB2/CB2 lie considerably lower/higher in energy and are therefore not seen in this figure. A band gap of 34 meV has opened up. The bands are very flat in the vicinity of K, and we therefore do not define a Dirac point.

Finally, we consider the band structure of the four kinds of SW defects in TBLG; these results, in the vicinity of K, are shown in Figs. cf. It is interesting to compare these figures with Fig. 4, which shows the band structure of TBLG in the absence of the SW defect. We see that the introduction of the SW defect has had a significant impact on the band structure of TBLG. Moreover, the four figures  cf show perceptible differences. This is because, as can be seen in Fig. 5, the registry between the atoms of the two layers is quite different for the four types of SW defects, as a result of which the interlayer coupling is quite different in the four cases.

Let us first consider Fig. c, which is the easiest of the four to understand in terms of how it arises from two twisted bilayers, one of which contains a SW defect. The Dirac points of the two layers are shifted from each other in k-space, even when we consider the “small” Brillouin zone that is commensurate with both layers [See the blue hexagon in Fig. 2, also shown in Fig. 7]. We wish to emphasize that this is different from the usual case of pristine TBLG with, e.g., the supercell, where the Dirac points corresponding to the upper and lower layers are shifted in extended k-space, but fold back to the K point when considering the “small” Brillouin zone corresponding to the commensurate supercell. For conceptual purposes, it is easy to think of each of these two Dirac points as arising from the Dirac points of the two individual layers, one unrotated and pristine, and the other rotated and containing a SW defect (note, however, that we will see further below that this picture is oversimplified, and needs to be further qualified). For TBLG with a SW defect, the Dirac point that arises primarily from the unrotated and undefective layer lies very close to the K point [see the position of the blue diamond that lies within the first BZ in Fig. 7(b)]. The second Dirac point, which arises primarily from the rotated defective layer, is shifted further away from K [see the position of the second blue diamond in Fig. 7(b)], as a result of the SW defect, as was seen, e.g., in Fig. a. Very small minigaps open up at the two Dirac points due to interactions between the two layers. The two Dirac crossing energies and are slightly shifted with respect to each other, with one lying slightly below the Fermi energy , and the other lying slightly above it. At energies further away from , the two Dirac cones intersect, and the resulting avoided crossings result in the opening up of band gaps, as a result of which VB1 has a skewed “M” shape, CB1 has a skewed “W” shape, and the bands VB2 and CB2 directly below and above these have conical shapes. The saddle points in VB1 and CB1 that form where the two cones fuse will give rise to van Hove singularities in the electronic density of states.

The same kind of interpretation applies to the other three types of SW defects (, and ) shown in Figs. 8(d)–(f), though in these cases the underlying double cone structure is more distorted since the gaps that open up are larger, because of a greater effect of interlayer coupling. We point out that one effect of the introduction of the SW defects is an increase in the density of states at the Fermi energy, since the Dirac points and have now been shifted slightly below and above .

As already noted, in Fig. 7(b) we have marked the positions in k-space of the Dirac points for the four kinds of SW defects in our TBLG. The Dirac points have been defined as the points in k-space corresponding to the maximum of VB1 and the minimum of CB1. Of these, one (arising primarily from the pristine layer) lies close to K, and one (arising primarily from the layer with the Stone-Wales defect) further away. The direction of the shift of the latter in k-space is determined by the orientation of the SW defect in real space.

The skewed shape of the band structure arising from double Dirac cones [see Fig. 8], where VB1 looks like a tilted “M”, and CB1 like a tilted “W”, is reminiscent of the effects of a transverse electric field on the band structure of twisted bilayer graphene.(4); (7) Note that even quite small electric fields (arising, e.g., due to a very minute charge transfer) can be expected to have a big effect on the positioning of the Dirac points with respect to , because of the very low electronic density of states in this region. We next check what the nature of the charge transfer (between the two graphene layers) and charge redistribution (in the vicinity of each graphene layer) is when two graphene sheets are brought together to form our TBLG systems. In Fig. 10 (a) we have plotted our results for , the planar integral of the change in electronic charge density of TBLG with a SW defect, referenced to the sum of the charge densities of an isolated pristine graphene layer and an isolated graphene layer with a SW defect. Note that the charge densities of the latter two systems are computed making use of their geometries in the combined system. For comparison, we also plot the same results for TBLG without a SW defect. One can clearly see that while the curve for TBLG with a SW defect (red solid line) is asymmetric about the midpoint, that for TBLG with no defect (black dashed line) is symmetric. By integrating outward from the midpoint, i.e., , we get the results shown in Fig. 10 (b). This shows clearly that upon bringing the two layers together, there is a depletion of electronic charge from the pristine layer, and an accumulation of charge in the defective layer; the net charge transferred is 0.0106 electrons. It is therefore quite tempting to apply a model where there is a uniform positive charge (depletion of electrons) on the pristine layer and uniform negative charge (accumulation of electrons) on the defective layer, and approximate its effects as those of an electric field due to a parallel plate capacitor. This would yield an electric field of strength NC. We see from Fig. 8(c) that in our case, shifts by 0.02 eV with respect to . These results appear to be consistent with those of a previous study of undefective TBLG in an electric field, where the authors found that an electric field of NC resulted in the Dirac points getting shifted with respect to by 0.1 eV. However, while these results might seem initially encouraging, the actual situation is not so simple, since in our case, the shift in energy of the two Dirac points is opposite in sign to that predicted by this simple model of charge transfer, whether interpreted in terms of the direction of the electric field, or whether interpreted (equivalently) in terms of how and should shift with respect to one another when electrons are removed or added.

A similar situation, in which the relative shift between and is, in some cases, in the direction opposite to that expected from a naive model of charge transfer, has been observed by previous authors who studied the electronic structure of graphene on metals.(66); (67); (68); (69) The authors of these previous studies attributed this to the presence of chemical interactions between the graphene and metal surface, which cannot be accounted for by the simple picture of charge transfer. For our systems, we find unmistakable signatures of the interaction between the two layers; it is clear that a scenario in which one Dirac cone is seen as arising entirely from the pristine layer, and the other as arising entirely from the defective layer, is not correct. As evidence of this, in Fig. 9, we have plotted the charge densities of four states very near the Dirac points, one each just above and below , and one each just above and below , for the case of the Stone-Wales defect in TBLG. In all four cases, it is obvious that the states have significant charge density on both graphene layers, even though the two states near have a greater charge density on the pristine (lower) layer, and those near have a greater charge density on the defective (upper) layer.

We have also computed the band-structure of the SW defect in TBLG using LDA. Upon comparing with the results obtained using DFT-D2, (see Fig. 8) we find almost no change in the vicinity of , though there are some slight differences at energies away from .

iv.3 Monovacancy Defect

A monovacancy defect can be created in a graphene sheet by removing one of the carbon atoms. In our example of TBLG, there are three inequivalent choices for the site at which this can be done; we term the resulting monovacancies , and (see Fig. 1). We consider one defect per supercell (corresponding to a defect density of 1.8%) for all these three types of monovacancies. For the case of the -monovacancy, we have also considered different supercells and , corresponding to defect densities of 7.2% and 0.8%, respectively. This allows us to examine the role of defect-defect interactions and strain relaxation in the defective systems.

To create these vacancies one needs to supply energy to the system; this is known as the monovacancy formation energy. The formation energy of a monovacancy is defined as:(70)

 ΔEvac=EvacSi−(N−1N)EnovacSi (14)

where is the total energy of the pristine system containing atoms in the supercell , and is the total energy of the supercell that contains a single monovacancy and is comprised of atoms. Note that the system with the monovacancy is found to possess a magnetic moment, and all these total energies are therefore obtained from spin-polarized DFT calculations.

Our results for the monovacancy formation energy for various single-layer and bilayer graphene systems are shown in Table 3. As in the case of the SW defect, we see that the different atomic arrangements in the different types of monovacancies are not significantly reflected in the energetics of defect formation: the values of for the three types of monovacancy for TBLG are very close to one another, and also to the values for SLG and AB-BLG.

Structural Properties

In all the three types of monovacancy considered for TBLG, both layers remain almost planar. The nearest-neighbor atoms of the vacancy site are displaced noticeably from their positions in the undefective structure. The removal of a carbon atom results in the under-coordination of its three nearest neighbor carbon atoms, resulting in the creation of a lone-pair of electrons over each carbon atom. The three-fold symmetry about the vacancy site is broken by a Jahn-Teller type distortion, resulting in an effective bonding between two of the three nearest neighbor carbon atoms. We find that two of the atoms are attracted toward each other and the distance between these atoms is reduced significantly, showing a signature of bonding. A similar phenomenon has earlier been reported for a monovacancy in a single layer of graphene (21) and for an AB-stacked bilayer of graphene.(71) We observe such a phenomenon in all the three vacancy positions: , and -vacancy [see Figs. 11, 11 and 11]. For the choice of defective supercell , this bond distance is found to be 2.05 Å, 2.03 Å, and 2.04 Å, for , and -vacancies, respectively. For the case of the -vacancy, we have also considered different supercells and , corresponding to a defect density of 7.2% and 0.8%, respectively. In these cases, this shorter bond distance is found to be 2.29 Å and 1.86 Å, respectively.

Simulated STM images

Next, we compute simulated STM images of the monovacancy defects in TBLG (See Fig. 11), making use of a bias voltage of 0.4 eV above the Fermi energy. The STM images of the three type of monovacancies (, and ) are basically indistinguishable, if the STM image is taken from the defective side of the bilayer. Interestingly, however if the STM image is taken from the undefective side of the bilayer, one can clearly distinguish between the three types of monovacancies, as their STM images show markedly different signatures; see the three lowest panels in Fig. 11.

Electronic Structure and Magnetic Properties

Creating a monovacancy in TBLG results in the appearance of a magnetic moment; its origin is similar to that in single layer graphene.(25); (23) In a sheet of pristine graphene, each carbon atom is bonded to three of its neighbors by three hybridized -bonds and one -bond, each bond sharing two electrons. When a monovacancy is created, the missing carbon atom takes away its share of four electrons, and there are three unsatisfied -electrons (one electron localized on each of the three nearest neighbor carbon atoms) and one -electron, in the vicinity of the monovacancy. Thus an electron in each of these states will try to maximize its spin, giving rise to a net magnetic moment of 4 . As already mentioned, of the three C atoms surrounding the vacancy, two atoms rebond, lowering the energy via a Jahn-Teller type distortion, and as a result the spins of the -electrons involved in this bonding orient antiparallel to each other, due to Pauli’s exclusion principle. Thus the monovacancy is expected to have a net magnetic moment of 2 , arising from one state and one state. This is why creating a monovacancy in a graphene sheet results in the appearance of a magnetic moment. However, as the -electrons have a somewhat itinerant character, and their bands cross the Fermi level, there is fractional occupation of this state so that the magnetic moment becomes less than 2 . In our case of a monovacacy in TBLG, in supercells of different sizes , and , we find a net magnetic moment of 1.81 , 1.25 and 1.85 per defect, respectively. These values are in good agreement with earlier reported values for monovacancy defects in SLG (magnetic moments of ranging from 1 to 2 depending on the defect concentration (21); (25)) and in AB-BLG (magnetic moment of 1.3 (71)). In Table 3, we have reported the values, for each of the cases considered by us in the supercell , for the net magnetic moment per defect defined as and the absolute magnetic moment per defect defined as . Here and are the up-spin and down-spin charge densities, respectively, and the integral is carried out over the supercell . All the three exchange-correlation functionals give for all systems. Interestingly however, the values obtained for vary with the functionals used, being for LDA and for PBE and DFT-D2. No significant changes are seen on going from SLG to AB-BLG or TBLG.

The effect of twist on the electronic structure of bilayer graphene with a monovacancy is manifested most clearly when one examines the band structure; these results are shown in Fig. 12. The presence of the vacancy in the twisted supercell results in an IBZ which is half of the entire BZ, similar to the case of the SW defect. As there is no shift in the position in k-space of the Dirac point, in the case of monovacancy defects, it suffices to plot the band structure along the conventionally used paths in the BZ; this is what we have done. We first compare the results obtained in the three types of monovacancy defects in TBLG [see Fig. 12]. We see that for a single type of monovacancy, the features of the four bands (two up-spin and two down-spin) that cross the Fermi level are largely similar in each third of the IBZ. Moreover, these four bands are also quite similar for the three different types of monovacancies (differences in the band structure become more apparent as one moves further away from the Fermi level). Therefore, for simplicity, when discussing the features of the band structure, and comparing it with the band structure of a monovacancy in SLG and pristine SLG, we will restrict ourselves to considering the band structure of the -monovacancy in TBLG, in one-third of its IBZ.

This is the comparison we carry out in Fig. 13. We see that the band structure of TBLG with a monovacancy [Fig. 13(c)] is clearly derived from the superposition of the band structures of a pristine SLG layer [Fig. 13(a)] and SLG with a monovacancy [Fig. 13(b)] with some important modifications. The bands derived from those of pristine SLG are shifted up in energy quite significantly, so that the Dirac crossing point now occurs eV above . In contrast, the bands arising from the layer with the monovacancy defect are shifted very slightly down in energy, there is also a slight change in their shape, which has a significant impact on the density of states because of the flat dispersion.

The shifting up/down of bands has three causes: (i) the fact that the layer with the monovacancy has four fewer valence electrons than the pristine layer (ii) a charge transfer between the layers, similar to that observed in the case of the Stone-Wales defect (iii) a chemical interaction between the layers. The effect of (i) alone can be artificially examined by separating the two layers by a very large distance of 6.5 Å, in that case, we find that the Dirac crossing energy arising from the pristine layer lies 0.32 eV above ; this shift can be attributed entirely to the requirement of equilibration of the Fermi energies of the two layers. However, when the spacing between the two layers is reduced to its equilibrium value of 3.3 Å, this upward shift is partially cancelled by a downward shift that arises due to the net combined effect of (ii) and (iii), resulting in the final upward shift of 0.25 eV observed. The initial shift of 0.32 eV can also be understood from computations of the work function . For SLG, we find that eV, whereas for SLG with a monovacancy, eV; both these values are computed making use of the same unit cell that is used for calculations of TBLG. When the two layers are brought together, with an interlayer separation of 6.5Å, we again obtain eV; this is because of the flat dispersion of the monovacancy states near and the resulting high value of the density of states, so that the states arising from the layer with the monovacancy do not shift perceptibly in energy, whereas the states arising from the pristine layer are shifted in energy by eV.

The charge transfer can be computed, as in the case of the Stone-Wales defect, by integrating the planar average of the charge redistribution. The planar average is plotted as a function of the coordinate in Fig. 10(c) and is seen to be asymmetric. On integrating outward from the midpoint of the two graphene layers, we find that the pristine layer has lost 0.02 electrons, and the layer with the monovacancy has gained 0.02 electrons [see Fig. 10(d)]. These values are in fairly good accordance with the values obtained from a Bader charge analysis, which yields a net charge on the pristine layer of 223.983 electrons, and on the defective layer of 220.016 electrons, i.e., 0.016 electrons have been transferred from the former to the latter.

The redistribution in charge upon bringing the two twisted layers (one pristine, and one with a monovacancy) together can be seen clearly upon plotting a charge density difference map (See Fig. 14). To obtain this map, we subtract the individual charge densities of the defective SLG and the pristine SLG (maintaining their geometry as that in the combined system), from the charge density of the defective bilayer. A region of charge depletion is developed below the vacancy site, slightly above the undefective layer. We have plotted this for the monovacancy, similar plots are obtained for the and monovacancies. The charge redistribution in the presence of a monovacancy is again a manifestation of the interlayer coupling in twisted bilayer graphene.

Once again, as in the case of TBLG with a Stone-Wales defect, if one were to consider merely the simple model of the net charge transfer from the undefective to the defective layer, one would have expected to shift up further on reducing the separation from 6.55 Å to 3.3 Å. The fact that, in contrast, it shifted downward, is a manifestation of the chemical interaction between the layers, as already discussed in the case of the SW defect.

The magnetic defect state that arises due to the creation of the vacancy in a single layer of graphene (SLG) can be seen in the band structure as a spin-polarized flat band about 0.5 eV below the Fermi energy , and can also be identified as a sharp peak in the density of states [See Figs. 13(b) and (c)]. This corresponds to the unsatisfied -bond on one of the atoms near the vacancy. Further, the monovacancy acts as a scattering center for the Dirac-like -electrons, thus resulting in the breakdown of the Dirac cone in the defective sheet, instead opening up a gap. The Fermi energy is lowered due to the removal of a carbon atom, and hence the Fermi level falls in the bands below this opened gap. The -band also becomes spin polarized, as can be seen in the energy range -0.25 eV to 0.4 eV. The high density of states near the Fermi energy for one of the spin channels leads to a high spin polarization of -electrons below .

We define the spin-polarized charge density as the difference between the majority and minority spin charge densities. In Fig. 15, we have plotted for states that fall within an energy window of width 0.05 eV below , for an monovacancy in TBLG. It is clearly evident that the spin polarization is not restricted to the immediate neighborhood of the monovacancy, but is widespread in spatial extent. This may possibly lead to a high degree of spin polarization of a current that is passed through such a system. However, in a situation of randomly oriented monovacancy defects, the net spin-polarization of the charge density will possibly depend on the sublattice details of a particular distribution of the defects, similar to that seen in the case of SLG.(23)

V Summary and Conclusions

In this study, we have considered two types of point defects, viz., Stone-Wales defects and monovacancy defects, that can arise or can be induced in twisted bilayer graphene. We have studied the effect of these defects on the structural, electronic and magnetic properties of twisted bilayer graphene. We have compared these results with the same defects in single layer graphene and AB-stacked bilayer graphene.

The formation energy of these defects is very similar to those of the corresponding defects in single layer graphene and in AB-stacked graphene; this is to be expected because the coupling between the two layers in twisted bilayer graphene is weak, and thus the effect on energetics of the second twisted layer is not particularly marked. Similarly, the defect formation energy does not change appreciably when the defect is induced on various symmetry-inequivalent sites on TBLG.

However, signatures of interlayer coupling are manifested clearly when examining the band structure of TBLG with defects. When we consider TBLG consisting of one pristine layer and one with a SW defect, one gets the formation of two Dirac cones, shifted slightly in energy and momentum. The shift in energy of the two Dirac cones arises due to the combined effect of a transfer of electrons from the pristine to the defective layer, and a chemical interaction between the two layers; this mimics the effect of placing the system in an electric field. The band structure depends on the site at which the SW defect is created – at some sites, the underlying double cone structure near the K points is distorted by the opening up of mini gaps – it might however be difficult to control the site at which a defect is created in an actual experimental situation.

We note here that the use of periodic boundary conditions in our plane wave calculations means that we actually have a periodic array of SW defects. Recent work has shown that it is possible to create multiple Dirac cones by artificially engineering a lateral periodic potential by creating a superlattice.(73)

Creating a monovacancy in one layer of TBLG shifts up the Dirac cone of the other layer by eV, due to a combination of the requirement of the equilbration of the Fermi energies, charge transfer and chemical interaction between the sheets. Both SW defects and monovacancies increase the density of states at ; this should be of relevance to transport properties. The monovacancy defect in TBLG leads to a magnetic defect state, very similar to the case of monovacancy defect in SLG and AB-BLG. A sharp magnetic defect state is the signature of the monovacancy, and arises from the dangling -bond on one of the nearest neighbor atoms of the vacant site. The defect has a magnetic moment of 1.2 – 2.0 , depending on the defect density. The scattering of the -states of the defective layer results in a highly spin-polarized density of states at the Fermi energy, both in the case of SLG and TBLG. The spin-polarized states in the vicinity of the Fermi energy are widespread in spatial extent, thus raising the possibility of a large spin-polarization of the current in the defective layer.

By comparing results obtained with the LDA, GGA (PBE) and DFT-D2 methods, we find that the PBE alone (i.e., without van der Waals interactions incorporated) is generally inadequate to describe bilayer graphene in general, and TBLG in particular. In general, both LDA and DFT-D2 give very similar results, with a few notable exceptions: DFT-D2 gives a value for the exfoliation energy that is in closer agreement with experiment, and the two methods, while giving similar results for the net magnetic moment , give different values for the absolute magnetic moment for TBLG with a monovacancy.

References

1. A. K. Geim and K. S. Novoselov, Nature Materials 6 183 (2007).
2. A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, Rev. Mod. Phys. 81, 109 (2009).
3. K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, M. I. Katsnelson, I. V. Grigorieva, S. V. Dubonos and A. A. Firsov, Nature (London) 438, 197 (2005).
4. J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H. Castro Neto, Phys. Rev. Lett. 99, 256802 (2007).
5. S. Shallcross, S. Sharma, and O. A. Pankratov, Phys. Rev. Lett. 101, 056803 (2008).
6. A. Luican, G. Li, A. Reina, J. Kong, R. R. Nair, K. S. Novoselov, A. K. Geim, and E. Y. Andrei, Phys. Rev. Lett. 106, 126802 (2011).
7. G. Li, A. Luican, J. M. B. Lopes dos Santos, A. H. Castro Neto, A. Reina, J. Kong and E. Y. Andrei, Nature Physics 6, 109 (2010).
8. P. Moon and M. Koshino, Phys. Rev. B 85, 195458 (2012).
9. L. Xian, Z. F. Wang, and M. Y. Chou, Nano Lett. 13, 5159 (2013).
10. J. Hass, F. Varchon, J. E. Millán-Otoya, M. Sprinkle, N. Sharma, W. A. de Heer, C. Berger, P. N. First, L. Magaud, and E. H. Conrad, Phys. Rev. Lett. 100, 125504 (2008).
11. L. Meng, Z. D. Chu, Y. Zhang, J. Y. Yang, R. F. Dou, J. C. Nie, and L. He, Phys. Rev. B 85, 235453 (2012).
12. Z. Y. Rong and P. Kuiper, Phys. Rev. B 48, 17427 (1993).
13. T. Ohta, J. T. Robinson, P. J. Feibelman, A. Bostwick, E. Rotenberg, and T. E. Beechem, Phys. Rev. Lett. 109, 186807 (2012).
14. A. Hashimoto, K. Suenaga, A. Gloter, K. Urita, and S. Iijima, Nature (London) 430, 870-873 (2004).
15. A. J. Stone and D. J. Wales, Chem. Phys. Lett. 128, 501 (1986).
16. F. Banhart, J. Kotakoski, and A. V. Krasheninnikov, ACS Nano 5, 26 (2011).
17. D. Zhan, J. Yan, L. Lai, Z. Ni, L. Liu, and Z. Shen, Adv. Mater. 24, 4055-4069 (2012).
18. L. Pisani, B. Montanari and N. M. Harrison, New Journal of Physics 10, 033002 (2008).
19. J-H. Chen, W. G. Cullen, E. D. Williams, and M. S. Fuhrer, Nature Physics 7, 535 (2011).
20. P. O. Lehtinen, A. S. Foster, A. Ayuela, A. Krasheninnikov, K. Nordlund, and R. M. Nieminen, Phys. Rev. Lett. 91, 017202 (2003).
21. P. O. Lehtinen, A. S. Foster, Y. Ma, A. V. Krasheninnikov, and R. M. Nieminen, Phys. Rev. Lett. 93, 187202 (2004).
22. O. V. Yazyev and L. Helm, Phys. Rev. B 75, 125408 (2007).
23. O. V. Yazyev, Rep. Prog. Phys. 73, 056501 (2010).
24. Y. Ma, P. O. Lehtinen, A. S. Foster and R. M. Nieminen, New Journal of Physics 6, 68 (2007).
25. B. R. K. Nanda, M. Sherafati, Z. S. Popović, and S. Satpathy, New Journal of Physics 14, 083004 (2012).
26. Y. Zhang, Y. Chen, K. Zhou, C. Liu, J. Zeng, H. Zhang and Y. Peng, Nanotechnology 20, 185504 (2009).
27. S. Y. Zhou, G. Gweon, A. V. Fedorov, P. N. First, W. A. De Heer, D. Lee, F Guinea, A. H. Castro Neto, and A. Lanzara, Nat. Mater. 6, 770 (2007).
28. X. Y. Peng and R. Ahuja, Nano Lett. 8, 4464 (2008).
29. Y. H. He, L. Wang , X. L. Chen, Z. F. Wu, W. Li, Y. Cai, and N. Wang, Appl. Phys. Lett. 99, 033109 (2011).
30. G. X. Liu, D. Teweldebrhan, and A. A. Balandin, IEEE Trans. Nanotechnol. 10, 865 (2011).
31. D. Teweldebrhan and A. A. Balandin, Appl. Phys. Lett. 94, 013101 (2009).
32. D. Teweldebrhan and A. A. Balandin, Appl. Phys. Lett. 95, 246102 (2009).
33. S. Mathew, T. K. Chan, D. Zhan, K. Gopinadhan, A. R. Barman, M. B. H. Breese, S. Dhar, Z. X. Shen, T. Venkatesan, and J. T. L. Thong, J. Appl. Phys. 110, 084309 (2011).
34. S. Mathew, T. K. Chan, D. Zhan, K. Gopinadhan, A. R. Barman, M. B. H. Breese, S. Dhar, Z. X. Shen, T. Venkatesan, and J. T. L. Thong, Carbon 49, 1720 (2011).
35. J. W. Bai, X. Zhong, S. Jiang, Y. Huang, and X. F. Duan, Nat. Nanotechnol. 5, 190 (2010).
36. M. M. Ugeda, I. Brihuega, F. Guinea, and J. M. Gomez-Rodriguez, Phys. Rev. Lett. 104, 096804 (2010).
37. S. Stankovich, D. A. Dikin, G. H. B. Dommett, K. M. Kohlhaas, E. J. Zimney, E. A. Stach, R. D. Piner, S. T. Nguyen, and R. S. Ruoff, Nature (London) 442, 282 (2006).
38. S. Ryu, M. Y. Han, J. Maultzsch, T. F. Heinz, P. Kim, M. L. Steigerwald, and L. E. Brus, Nano Lett. 8, 4597 (2008).
39. D. C. Elias, R. R. Nair, T. M. G. Mohiuddin, S. V. Morozov, P. Blake, M. P. Halsall, A. C. Ferrari, D. W. Boukhvalov, M. I. Katsnelson, A. K. Geim, and K. S. Novoselov, Science 323, 610 (2009).
40. Z. Q. Luo, J. Z. Shang, S. H. Lim, D. H. Li, Q. H. Xiong, Z. X. Shen, J. Y. Lin, and T. Yu, Appl. Phys. Lett. 97, 233111 (2010).
41. Z. H. Ni, L. A. Ponomarenko, R. R. Nair, R. Yang, S. Anissimova, I. V. Grigorieva, F. Schedin, P. Blake, Z. X. Shen, E. H. Hill, K. S. Novoselov, and A. K. Geim, Nano Lett. 10, 3868 (2010).
42. R. R. Nair, W. C. Ren, R. Jalil, I. Riaz, V. G. Kravets, L. Britnell, P. Blake, F. Schedin, A. S. Mayorov, S. J. Yuan, M. I. Katsnelson, H. M. Cheng, W. Strupinski, L. G. Bulusheva, A. V. Okotrub, I. V. Grigorieva, A. N. Grigorenko, K. S. Novoselov, and A. K. Geim, Small 6, 2877 (2010).
43. K. J. Jeon, Z. Lee, E. Pollak, L. Moreschini, A. Bostwick, C. M. Park, R. Mendelsberg, V. Radmilovic, R. Kostecki, T. J. Richardson, and E. Rotenberg, ACS Nano 5, 1042 (2011).
44. P. Giannozzi, et al., J. Phys.: Condens. Matter 21, 395502 (2009).
45. D. Vanderbilt, Phys. Rev. B 41, 7892 (1990).
46. J. P. Perdew and A. Zunger, Phys. Rev. B 23, 5048 (1981).
47. J. P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 77,3865-3868 (1996).
48. S. Grimme, J. Comput. Chem. 25, 1463 (2004).
49. S. Grimme, J. Comput. Chem. 27, 1787 (2006).
50. N. Marzari, D. Vanderbilt, A. De Vita, and M. C. Payne, Phys. Rev. Lett. 82, 3296 (1999).
51. R. P. Feynman, Phys. Rev. 56 340 (1939).
52. J. Tersoff and D. R. Hamann, Phys. Rev. Lett. 50, 1998 (1983).
53. S. Shallcross, S. Sharma, E. Kandelaki, and O. A. Pankratov, Phys. Rev. B 81, 165105 (2010).
54. J. D. Bernal, Proc. R. Soc. London Ser. A 106, 749 (1924).
55. L. A. Girifalco and R. A. Ladd, J. Chem. Phys. 25, 693 (1956).
56. L. X. Benedict, N. G. Chopra, M. L. Cohen, A. Zettl, S. G. Louie, and V. H. Crespi, Chem. Phys. Lett. 286, 490 (1998).
57. Z. Liu, J. Z. Liu, Y. Cheng, Z. Li, L. Wang, and Q. Zheng, Phys. Rev. B 85, 205418 (2012).
58. S. Latil, V. Meunier, and L. Henrard, Phys. Rev. B 76 201402(R) (2007).
59. Y. Mao and J. Zhong, Nanotechnology 19, 205708 (2008).
60. G. M. Rutter, S. Jung, N. N. Klimov, D. B. Newell, N. B. Zhitenev, and J. A. Stroscio, Nature Physics 7, 649 (2011).
61. J. Ma, D. Alfè, A. Michaelides, and E. Wang, Phys. Rev. B 80, 033407 (2009).
62. L. G. Zhou, and S. Shi, Appl. Phys. Lett. 83, 1222-1224 (2007).
63. V. N. Popov, L. Henrard, and P. Lambin, Carbon 47, 2448 (2009).
64. S. N. Shirodkar and U. V. Waghmare, Phys. Rev. B 86, 165401 (2012).
65. S. Pisana, M. Lazzeri, C. Casiraghi, K. S. Novoselov, A. K. Geim, A. C. Ferrari, and F. Mauri, Nat. Mater. 6, 198 (2007).
66. G. Giovannetti, P. A. Khomyakov, G. Brocks, V. M. Karpan, J. van den Brink, and P. J. Kelly, Phys. Rev. Lett. 101, 026803 (2008).
67. P. A. Khomyakov, G. Giovannetti, P. C. Rusu, G. Brocks, J. van den Brink, and P. J. Kelly, Phys. Rev. B 79, 195425 (2009).
68. C. Gong, G. Lee, B. Shan, E. M. Vogel, R. M. Wallace, and K. Chao, J. Appl. Phys. 108, 123711 (2010).
69. K. Pi, K. M. McCreary, W. Bao, W. Han, Y. F. Chiang, Y. Li, S.-W. Tsai, C. N. Lau, and R. K. Kawakami, Phys. Rev. B 80, 075406 (2009).
70. R. Singh and P. Kroll, J. Phys.: Condens. Matter 21, 196002 (2009).
71. S. Choi, B. Jeong, S. Kim and G. Kim, J. Phys. Cond. Matt. 20, 235220 (2008).
72. P. A. Thrower and R. M. Mayer, Phys. Stat. Sol. A 47 11 (1978).
73. S. Dubey, V. Singh, A. K. Bhat, P. Parikh, S. Grover, R. Sensarma, V. Tripathi, K. Sengupta, and M. M. Deshmukh, Nano Lett. 13, 3990 (2013).
180858