week ending 7 MARCH 2014
PHYSICAL REVIEW LETTERS
PRL 112, 096805 (2014)
Strain-Induced Pseudomagnetic Fields in Twisted Graphene Nanoribbons 1
Dong-Bo Zhang,1,4 Gotthard Seifert,2 and Kai Chang3,*
Department of Chemical Engineering and Materials Science, University of Minnesota, 421 Washington Avenue SE, Minneapolis, Minnesota 55455, USA 2 Physikalische Chemie, Technische Universität Dresden, D-01062 Dresden, Germany 3 SKLSM, Institute of Semiconductors, Chinese Academy of Sciences, P.O. Box 912, Beijing 100083, China 4 Beijing Computational Science Research Center, Beijing 100084, China (Received 29 May 2013; published 6 March 2014) We present, for the first time, an atomic-level and quantitative study of a strain-induced pseudomagnetic field in graphene nanoribbons with widths of hundreds of nanometers. We show that twisting strongly affects the band structures of graphene nanoribbons with arbitrary chirality and generates well-defined pseudo-Landau levels, which mimics the quantization of massive Dirac fermions in a magnetic field up to 160 T. Electrons are localized either at ribbon edges forming the edge current or at the ribbon center forming the snake orbit current, both being valley polarized. Our result paves the way for the design of new graphene-based nanoelectronics. DOI: 10.1103/PhysRevLett.112.096805
PACS numbers: 73.22.-f, 73.23.-b, 73.40.Gk, 85.30.De
Owing to the linear dispersion and the chiral property of electrons at K and K 0 valleys where the charge carriers present the characteristic of massless Dirac fermions [1,2], graphene possesses a wealth of new physics such as the room temperature quantum Hall effect [3], tunable electron mobility [4], valley polarization [5,6], and furthermore the controllable electron-structure coupling [7]. Modeling of structural deformation with a gauge field [8–15] reveals that strong pseudomagnetic fields can be induced with a designed strain pattern [16–19] as demonstrated by experiment [20]. For these studies, it was found that the electronic states of deformed graphene undergo a Landau quantization, mimicking the effect of a real magnetic field [1,3,21–23]. Such insight renews the interest of strain engineering of graphene. Indeed, designing strain aligned along specific crystallographic directions to generate strong gauge fields acting as a homogeneous pseudomagnetic field has been developing as a paradigm for tuning electronic quantum states [24]. In graphene, the pseudomagnetic field can approach hundreds of tesla [20], strong enough to result in an insulating bulk and a pair of counter-propagating edge currents arising from K and K 0 valleys, in analogue to the helical edge states in a two-dimensional topological insulator [25]. Mechanically, the graphene monolayer can sustain elastic strain up to 25% [26,27], allowing for diverse strain designs to functionalize graphene [28]. In fact, there have already been several efforts devoted to generating a pseudomagnetic field employing different strain forms such as bubbling [17] and bending [18,19]. Experimental realization of these strains relies on the patterned substrate, which poses a considerable challenge. In this Letter, we suggest a promising approach to create a pseudomagnetic field which induces a pure valley edge current and snake orbit in a twisted graphene nanoribbon (GNR). Such a twisting 0031-9007=14=112(9)=096805(5)
technique was successfully used in carbon nanotubes before [29,30]. The pseudomagnetic field generated in a twisted GNR can be derived by the gauge representation of electronic states [17–19,31,32]. Figure 1(a) displays a twisted GNR and its horizontal projection. It has a nearly quadratic spatial distribution of the normal strain ϵyy along the width dimension, x, Fig. 1(b). Therefore, the corresponding pseudomagnetic field for electronic states at the K valley of a GNR with chirality χ at twist rate γ is (see the Supplemental Material [33] for detailed information), BK ¼ C0 ððcos2 χ − 3sin2 χÞ cos χÞγ 2 x.
(1)
The coefficient C0 ¼ cβt0 ϕ0 a−1 0 , where t0 and a0 are the nearest neighbor hopping and the equilibrium C-C bond length of graphene, respectively. c and β [17] are both dimensionless parameters. ϕ0 is the magnetic flux quantum. It is easy to see that BK is anisotropic [34], being strongest at χ ¼ 0° (zigzag GNR), gradually weakening at 0° < χ < 30° (chiral GNR), and completely vanishing at χ ¼ 30° (armchair GNR), Fig. 1(c). On the other side, BK also varies along the width, x and changes sign at the GNR center (x ¼ 0). The pseudomagnetic field picture is instructive in elucidating the physical origin of the quantum Hall effect (QHE) and other electromagnetic properties of deformed graphene [35]. Notice that Eq. (1) is based on the nearest neighbor approximation π-orbital tight-binding model for weakly deformed graphene systems. However, generating a strong homogeneous pseudomagnetic field over a wide area usually involves a complex strain pattern and high strain level. Consequently, the dependence of π-electron hoppings between the nearest neighbors on the strain tensor ϵij ði; j ¼ x; yÞ is no longer linear, where the coupling of the electronic structure with
096805-1
© 2014 American Physical Society
PRL 112, 096805 (2014)
PHYSICAL REVIEW LETTERS
FIG. 1 (color online). (a) A twisted GNR (upper panel) and its horizontal projection (lower panel). The coordinate system is defined as x along the GNR width dimension and y along the GNR axis. Twisting is around the GNR axis. The twist-induced pseudomagnetic field of K valley-specified electronic states, BK directs reversely on the two sides of GNR axis. The pseudomagnetic field of K 0 valley electronic states BK0 ¼ −BK by time reversal symmetry. The edge currents JK and JK0 arising from K valley and K 0 valley, respectively, are counterpropagating. The current SK and SK0 associated with the snake orbit are also counterpropagating but in a serpentine manner around the GNR axis. (b) Strain distribution ϵyy along the GNR width, x at twist rate γ ¼ 0.61°=nm and (c) the corresponding pseudomagnetic field BK for 170 nm wide GNRs with different chiralities χ according to Eq. (1) with a0 ¼ 1.42 Å, β ¼ 2, c ¼ 1, and t0 ¼ 2.7 eV. (d) A near-field scanning optical microscope setup including contacts I (injector) and C (collector) for the detection of localized edge states (solid curve) and snake states (dashed curve).
elastic strain is not well defined [36,37]. Furthermore, not only nearest neighbor interactions determine the electronic structure, but also far-reaching interactions play a nonnegligible role [38]. These issues call for the investigation of the electronic states beyond the scope of the low-energy k · p Hamiltonian [2]. Atomistic quantum mechanical simulations represent a reliable alternative. Because the graphene systems used in QHE study are often of the dimension greater than 102 nm [3,17,21], having a few million atoms, full first principle and systematic quantum mechanical (QM) calculations formulated with periodic boundary conditions are impossible given the current computational capacity. Here, we showcase a study of strain-induced Landau quantization in graphene by realistic QM atomistic simulation of twisted GNRs with a width greater than 102 nm. The simulation is carried out with a density-functional theory based tight-binding (DFTB) theory [39,40]. Such a QM study is possible due to a symmetry-adapted scheme for the electronic systems [41]. Through this scheme, one can utilize the helical symmetry created by the twist to enable the realistic QM modeling of twisted GNR on a
week ending 7 MARCH 2014
basic repeating motif which contains the same N atoms in the primitive cell of the flat GNR. For example, for a zigzag GNR with a width of 170 nm, N ¼ 1608 (see Supplemental Material [33]). We have carried out DFTB calculations on a series of twisted GNRs with two types (zigzag and armchair) of edges. We optimize GNR atomistic geometries via a conjugate gradient energy minimization to the repeating motif. In our calculation, the greatest width of the GNR is about 176 nm comparable with the size of the sample used in QHE experiments [4,21]. Our result reveals a quite different picture of the Landau quantization in twisted GNRs from the one predicted by theory [17–19,31,32] as shown in Eq. (1). That is, the twist induces a well-defined electronic Landau spectrum not only in zigzag GNR but also in armchair GNR and, thus, forms localized edge states. Meanwhile, a snake orbit presents at the GNR center due to the opposite directing pseudomagnetic field on the two sides of GNR. Therefore, if any bias electrical potential is applied along the GNR axis, one shall be able to observe valley-specified edge currents, JK (JK0 ) and snake orbit currents, SK (SK0 ), schematically shown in Fig. 1(a). Figure 2 displays the band structure and density of states (DOS) of a twisted zigzag GNR of width 170 nm at (a) γ ¼ 0°=nm, and (b) γ ¼ 0.61°=nm. A very wide GNR mimics closely the band structure of bulk graphene, showing neatly linear dispersion for the valley-specified electronic states, Fig. 2(a). The extra bands associated with the zigzag edges, being dispersion free in the simple tight binding modeling with the nearest neighbor approximation, are now dispersive without the particle-hole symmetry. Mostly, this owes to the interaction between the next nearest neighbor carbon atoms [38]. Consequently, their contribution to the DOS is off the Dirac point, as indicated by the down arrow. Under twisting, we witness the onset of the flattening of valley-specified electronic bands and the DOS curve starts to fluctuate (see the online supplementary material). Thoroughly flatting of valley electronic bands starts at γ ¼ 0.61°=nm, producing sharp and separated peaks in the DOS curve, Fig. 2(b). These peaks correspond to the discrete Landau levels (LLs). It should be noted that the zero-energy LL at the Dirac point is formed from valley-specified electronic states, not from the electronic states arising from zigzag edges which, located between the zero and the first negative LLs in the DOS curve as indicated by down arrow, are not sensitive to twist. This is because they are mostly confined at very narrow regions around the GNR margins (see the Supplemental Material [33]). The twist-induced pseudomagnetic field only perturbs them weakly. Further twisting the GNR would not bring significant change in the DOS (see the Supplemental Material [33]). This means that there is a critical twist γ c , or equivalently, a saturated effective pseudomagnetic field, Beff . For this exemplified GNR, γ c ¼ 0.61°=nm. Thus, for γ > γ c , the obtained Landau spectrum stays stationary. More insight can be obtained from the fact that the effect of pseudomagnetic field on electronic states is in similar
096805-2
PRL 112, 096805 (2014)
PHYSICAL REVIEW LETTERS
FIG. 2 (color online). Electronic band structures (upper panel) and density of states (lower panel) of a 170 nm-wide zigzag GNR at twist rate (a) γ ¼ 0°=nm and (b) γ ¼ 0.61°=nm. The down arrows indicate the extra electronic states arising from zigzag edges. In (b), the solid dot (upper panel) indicates the snake orbit states, and the vertical dashed lines (lower panel) indicate the Landau levels with n ¼ 0, 1, 2, 3, 4 according to Eq. (2) with Beff ¼ 85 T. (c) The first Landau level energy, E1 , and (d) the effective magnetic field, Beff , versus GNR width W 0 .
week ending 7 MARCH 2014
A careful analysis of the experimental conditions reveals that the graphene bubbling is of a relatively small size of about only 10 nm. Over such a small size scale, it is possible to achieve a large strain gradient without requiring extremely high strain. This is consistent with our observation of larger Beff in narrower GNRs. Next we turn to the study of an armchair GNR. The theoretical prediction [17–19,31,32], as shown in Eq. (1) for twisted armchair GNRs, suggests that there would not be any Landau quantization since BK ¼ 0. Here, our result represents an exception of theory. For the flat case (γ ¼ 0°=nm), an armchair GNR displays the typical linear dispersion of electronic band structure and smooth DOS, Fig. 3(a). Without the presence of a zigzag edge, there are not any extra bands. Under twisting, the electronic states undergo the Landau quantization. Figure 3(b) shows the band structure and DOS at twist rate γ c ¼ 0.66°=nm. Again, the distribution of major DOS peaks agrees well with Eq. (2) with Beff ¼ 72 T. The nonvanishing pseudomagnetic field mainly comes from the distortion of the hexagonal lattice. Note that by time reversal symmetry, the twist-induced Landau states at K 0 valley is degenerate with those at K valley. Therefore, the corresponding pseudomagnetic field is simply −Beff . The large split of zero LL is related to the chirality property. There is also a small splitting of zero LL in twisted zigzag GNRs due to symmetry breaking of the A-B lattice. See the Supplemental Material [33] for an analysis of the symmetry breaking of the A-B lattice. A GNR mimics a Hall bar [3,21]. Because of the finite size in the width dimension, the valley-specified electronic states are pushed against the GNR margins given the strong pseudomagnetic field, Beff , forming the so-called edge
with that of a real magnetic field applied perpendicular to the graphene sheet. The formed Landau spectrum is [3,21] En ¼ sgnðnÞvF
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2eℏBeff jnj;
(2)
where the integer n indexes LLs. Indeed, with Beff ¼ 85 T, the major DOS peaks fall onto the vertical dashed lines representing the first few LLs defined by Eq. (2), Fig. 2(b). For a systematic study, we have carried out DFTB calculations on a series of GNRs with different widths. For each GNR, we determine γ c and Beff according to Eq. (2). For the evolution of the DOS with width, see the Supplemental Material [33]. We characterize the size dependence of the lowest LL and Beff by a simple fitting of the atomistic data, as E1 ∝ 1=W 0 , Fig. 2(c), and Beff ∝ 1=W 20 , Fig. 2(d), respectively. In the previous theoretical studies [17–19], the strain induced pseudomagnetic field is of tens of tesla. We suggest that this is mainly because the strain level is relatively low in these theoretical studies with the maximum strain ∼10%. Here, Beff is much larger, especially for narrow GNRs, Fig. 2(d), close to the experimental observation [20] where the strain (bubbling) induced pseudomagnetic field can be as large as 300 T.
FIG. 3 (color online). Electronic band structures (upper panel) and density of states (lower panel) of a 176 nm-wide armchair GNR at twist rate (a) γ ¼ 0°=nm and (b) γ ¼ 0.66°=nm. In (b), the solid dot (upper panel) indicates the snake states, and the vertical dashed lines (lower panel) indicate the LLs n ¼ 1, 2, 3, 4 according to Eq. (2) with Beff ¼ 136 T. The down arrows indicate the zero LL.
096805-3
PRL 112, 096805 (2014)
PHYSICAL REVIEW LETTERS
week ending 7 MARCH 2014
FIG. 5 (color online). Wave function spatial distributions (jΨj2 ) along armchair GNR width x (a) for the first four LLs n ¼ 0ðzero LLÞ, 1, 2, 3 as shown in Fig. 3(b) at κ ¼ π and (b) for the snake orbit whose location in the energy spectrum is indicated in Fig. 3(b) at κ ¼ 0.2π. In (a), the electronic states represented by the red dashed line and blue solid line are degenerate. The vertical dashed lines indicate the maxima of the wave functions for the different LLs.
FIG. 4 (color online). Wave function spatial distributions (jΨj2 ) along the zigzag GNR width x (a) for the first four LLs n ¼ 0, 1, 2, 3 as shown in Fig. 2(b) at κ ¼ 0 and (b) for the snake orbit whose location in the energy spectrum is indicated in Fig. 2(b) at κ ¼ 0.7π. In (a), the electronic states represented by the red dashed line and blue solid line are degenerate. The vertical dashed lines indicate the maxima of the wave functions for the different LLs.
states [42]. Figure 4(a) and Fig. 5(a) display the spatial distribution of wave functions of edge states along the GRN width for the studied zigzag GNR and armchair GNR, respectively. Indeed, these edge states are localized in the narrow area at the GNR margins and have different distributions for different LLs. Thus, for lower order (higher order) LLs, the valley edge currents transport via narrower (wider) atomic channels. The presence of a zigzag edge as in zigzag and chiral GNRs also brings extra edge states (see Supplemental Material [33]) which, however, are susceptible to edge disorder such as defect and reconstruction. The effects, such as elastic backscattering, are the consequences. The edge states completely insusceptible of edge disorder can be obtained in twisted armchair GNRs. This is because these states absolutely originate from valley-specified electronic states since there is no longer a zigzag edge presented, thus forming pure valley edge currents. On the other side, the oppositely directing pseudomagnetic field, Beff , on the two sides of the GNR results in a peculiar state—a snake orbit [43–45]— which is otherwise nonexistent in bent GNRs [18,19]. Spatially, it is located at the GNR center having a different momentum from edge states. See Fig. 4(b) and Fig. 5(b) for
the snake orbits in zigzag and armchair GNRs, respectively. To visualize such edge states and the snake state, a detector can be proposed utilizing state-of-art near-field scanning microwave microscopy [46,47], as schematically shown in Fig. 1(d). The twist strain can be applied in a similar way of the carbon nanotube pedal experiments [29,30]. The GNR is placed on a substrate of the insulating polymer. At the emergence of the edge current (snake current), which can be controlled by tuning the momentum of valley polarized electron beam, the highlight of the microwave image will occur at the margins (at the center) of the GNR. To summarize, we accomplished the realistic QM atomistic study of strain-induced electronic Landau quantization in GNRs under twisting, a fundamental type of deformation not considered before. Although the low energy theory [17–19,31,32] predicts that twisting is not an ideal candidate to induce Landau quantization since the strain induced pseudomagnetic field is essentially inhomogeneous (In particular it is zero for armchair GNR) as shown in Fig. 1(c), our calculation reveals that a well-defined Landau spectrum can be generated for both types (zigzag and armchair) of GNRs by twisting, demonstrating the importance of realistic atomistic calculation. We then immediately conclude that Landau quantization can also be achieved in other chiral GNRs. This revelation alleviates greatly the concern in experiment of controlling GNR’s orientation (chirality). The effective pseudomagnetic field orienting oppositely on the two sides of GNR, approaches hundreds of tesla. The electronic states, localized either at the GNRs edges, undergo Landau quantization, forming a valley-specified edge current, or at the GNR center, forming a valleyspecified snake orbit current. Note that in armchair GNR, we have pure valley edge currents. A doping effect [44] can
096805-4
PRL 112, 096805 (2014)
PHYSICAL REVIEW LETTERS
also be induced by twisting. Altogether, our results demonstrate a viable manner to control the electronic levels and transport properties of graphene at nanoscale, shedding new light on the design of novel electronics devices by strain engineering in graphene and other two-dimensional atomic crystals, e.g., silicene, MoS2 , etc. D.-B. Z. thanks Traian Dumitrica for discussion at the early stage of this work. G. S. acknowledges the support of the ERC project INTIF 226639. K. C. is supported by NSFC No. 10934007 and the National Basic Research Program of China (973 Program) under Grant No. 2011CB922204. Synergetic Innovation Center of Quantum Information and Quantum Physics at USTC, and CAS Center for excellence in quantum information and quantum physics.
*
[1]
[2] [3]
[4]
[5] [6] [7] [8] [9]
[10] [11] [12] [13] [14] [15] [16] [17] [18] [19]
Corresponding author. [email protected] 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). A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, Rev. Mod. Phys. 81, 109 (2009). K. S. Novoselov, Z. Jiang, Y. Zhang, S. V. Morozov, H. L. Stormer, U. Zeitler, J. C. Maan, G. S. Boebinger, P. Kim, and A. K. Geim, Science 315, 1379 (2007). K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A. Firsov, Science 306, 666 (2004). D. Xiao, W. Yao, and Q. Niu, Phys. Rev. Lett. 99, 236809 (2007). A. Rycerz, J. Tworzydo, and C. W. J. Beenakker, Nat. Phys. 3, 172 (2007). V. M. Pereira and A. H. Castro Neto, Phys. Rev. Lett. 103, 046801 (2009). M. A. H. Vozmediano, M. I. Katsnelson, and F. Guinea, Phys. Rep. 496, 109 (2010). S. V. Morozov, K. S. Novoselov, M. I. Katsnelson, F. Schedin, L. A. Ponomarenko, D. Jiang, and A. K. Geim, Phys. Rev. Lett. 97, 016801 (2006). M. M. Fogler, F. Guinea, and M. I. Katsnelson, Phys. Rev. Lett. 101, 226804 (2008). K. Sasaki, Y. Kawazoe, and R. Saito, Prog. Theor. Phys. 113, 463 (2005). A. F. Morpurgo and F. Guinea, Phys. Rev. Lett. 97, 196804 (2006). F. de Juan, A. Cortijo, and M. A. H. Vozmediano, Phys. Rev. B 76, 165409 (2007). F. Guinea, M. I. Katsnelson, and M. A. H. Vozmediano, Phys. Rev. B 77, 075422 (2008). E.-A. Kim and A. H. Castro Neto, Europhys. Lett. 84, 57007 (2008). F. Guinea, B. Horovitz, and P. Le Doussal, Phys. Rev. B 77, 205421 (2008). F. Guinea, M. I. Katsnelson, and A. K. Geim, Nat. Phys. 6, 30 (2009). F. Guinea, A. K. Geim, M. I. Katsnelson, and K. S. Novoselov, Phys. Rev. B 81, 035408 (2010). T. Low and F. Guinea, Nano Lett. 10, 3551 (2010).
week ending 7 MARCH 2014
[20] N. Levy, S. A. Burke, K. L. Meaker, M. Panlasigui, A. Zettl, F. Guinea, A. H. Castro Neto, and M. F. Crommie, Science 329, 544 (2010). [21] Y. Zhang, Y. -W. Tan, H. L. Stormer, and P. Kim, Nature (London) 438, 201 (2005). [22] D. A. Abanin, K. S. Novoselov, U. Zeitler, P. A. Lee, A. K. Geim, and L. S. Levitov, Phys. Rev. Lett. 98, 196806 (2007). [23] K. I. Bolotin, F. Ghahari, M. D. Shulman, H. L. Stormer, and P. Kim, Nature (London) 462, 196 (2009). [24] M. C. Rechtsman, J. M. Zeuner, A. Tünnermann, S. Nolte, M. Segev, and A. Szameit, Nat. Photonics 7, 153 (2012). [25] C. L. Kane and E. J. Mele, Phys. Rev. Lett. 95, 226801 (2005). [26] C. Lee, X. Wei, J. W. Kysar, and J. Hone, Science 321, 385 (2008). [27] K. S. Kim, Y. Zhao, H. Jang, S. Y. Lee, J. M. Kim, K. S. Kim, J. -H. Ahn, P. Kim, J. -Y. Choi, and B. H. Hong, Nature (London) 457, 706 (2009). [28] Z. H. Wu, F. Zhai, F. M. Peeters, H. Q. Xu, and K. Chang, Phys. Rev. Lett. 106, 176802 (2011). [29] T. Cohen-Karni, L. Segev, O. Srur-Lavi, S. R. Cohen, and E. Joselevich, Nat. Nanotechnol. 1, 36 (2006). [30] A. M. Fennimore, T. D. Yuzvinsky, W. -Q. Han, M. S. Fuhrer, J. Cumings, and A. Zettl, Nature (London) 424, 408 (2003). [31] J. L. Mañes, Phys. Rev. B 76, 045430 (2007). [32] H. Suzuura and T. Ando, Phys. Rev. B 65, 235412 (2002). [33] See Supplemental Material at http://link.aps.org/ supplemental/10.1103/PhysRevLett.112.096805 for the calculation method of twisted GNRs and detailed information about the spatial distributions of edge and snake orbit states. [34] T. O. Wehling, A. V. Balatsky, A. M. Tsvelik, M. I. Katsnelson, and A. I. Lichtenstein, Europhys. Lett. 84, 17003 (2008). [35] F. de Juan, A. Cortijo, M. A. H. Vozmediano, and A. Cano, Nat. Phys. 7, 810 (2011). [36] A. L. Kitt, V. M. Pereira, A. K. Swan, and B. B. Goldberg, Phys. Rev. B 85, 115432 (2012). [37] Y. Chang, T. Albash, and S. Haas, Phys. Rev. B 86, 125402 (2012). [38] N. M. R. Peres, F. Guinea, and A. H. Castro Neto, Phys. Rev. B 73, 125411 (2006). [39] D. Porezag, Th. Frauenheim, Th. Kohler, G. Seifert, and R. Kaschner, Phys. Rev. B 51, 12 947 (1995). [40] R. Ruralia and E. Hernández, Comput. Mater. Sci. 28, 85 (2003). [41] This special method was pioneered by several authors: C. T. White, D. H. Robertson, and J. W. Mintmire, Phys. Rev. B 47, 5485 (1993); V. N. Popov, New J. Phys. 6, 17 (2004); P. B. Allen, Nano Lett. 7, 1220 (2007). Also see Supplemental Material [33] for more details. [42] D. A. Abanin, K. S. Novoselov, U. Zeitler, P. A. Lee, A. K. Geim, and L. S. Levitov, Phys. Rev. Lett. 98, 196806 (2007). [43] J. E. Müller, Phys. Rev. Lett. 68, 385 (1992). [44] T. Low, F. Guinea, and M. I. Katsnelson, Phys. Rev. B 83, 195436 (2011). [45] J. R. Williams and C. M. Marcus, Phys. Rev. Lett. 107, 046602 (2011). [46] K. Lai, M. Nakamura, W. Kundhikanjana, M. Kawasaki, Y. Tokura, M. A. Kelly, and Z. -X. Shen, Science 329, 190 (2010). [47] K. Lai, W. Kundhikanjana, M. A. Kelly, Z. -X. Shen, J. Shabani, and M. Shayegan, Phys. Rev. Lett. 107, 176809 (2011).
096805-5
PHYSICAL REVIEW LETTERS
PRL 112, 096805 (2014)
Strain-Induced Pseudomagnetic Fields in Twisted Graphene Nanoribbons 1
Dong-Bo Zhang,1,4 Gotthard Seifert,2 and Kai Chang3,*
Department of Chemical Engineering and Materials Science, University of Minnesota, 421 Washington Avenue SE, Minneapolis, Minnesota 55455, USA 2 Physikalische Chemie, Technische Universität Dresden, D-01062 Dresden, Germany 3 SKLSM, Institute of Semiconductors, Chinese Academy of Sciences, P.O. Box 912, Beijing 100083, China 4 Beijing Computational Science Research Center, Beijing 100084, China (Received 29 May 2013; published 6 March 2014) We present, for the first time, an atomic-level and quantitative study of a strain-induced pseudomagnetic field in graphene nanoribbons with widths of hundreds of nanometers. We show that twisting strongly affects the band structures of graphene nanoribbons with arbitrary chirality and generates well-defined pseudo-Landau levels, which mimics the quantization of massive Dirac fermions in a magnetic field up to 160 T. Electrons are localized either at ribbon edges forming the edge current or at the ribbon center forming the snake orbit current, both being valley polarized. Our result paves the way for the design of new graphene-based nanoelectronics. DOI: 10.1103/PhysRevLett.112.096805
PACS numbers: 73.22.-f, 73.23.-b, 73.40.Gk, 85.30.De
Owing to the linear dispersion and the chiral property of electrons at K and K 0 valleys where the charge carriers present the characteristic of massless Dirac fermions [1,2], graphene possesses a wealth of new physics such as the room temperature quantum Hall effect [3], tunable electron mobility [4], valley polarization [5,6], and furthermore the controllable electron-structure coupling [7]. Modeling of structural deformation with a gauge field [8–15] reveals that strong pseudomagnetic fields can be induced with a designed strain pattern [16–19] as demonstrated by experiment [20]. For these studies, it was found that the electronic states of deformed graphene undergo a Landau quantization, mimicking the effect of a real magnetic field [1,3,21–23]. Such insight renews the interest of strain engineering of graphene. Indeed, designing strain aligned along specific crystallographic directions to generate strong gauge fields acting as a homogeneous pseudomagnetic field has been developing as a paradigm for tuning electronic quantum states [24]. In graphene, the pseudomagnetic field can approach hundreds of tesla [20], strong enough to result in an insulating bulk and a pair of counter-propagating edge currents arising from K and K 0 valleys, in analogue to the helical edge states in a two-dimensional topological insulator [25]. Mechanically, the graphene monolayer can sustain elastic strain up to 25% [26,27], allowing for diverse strain designs to functionalize graphene [28]. In fact, there have already been several efforts devoted to generating a pseudomagnetic field employing different strain forms such as bubbling [17] and bending [18,19]. Experimental realization of these strains relies on the patterned substrate, which poses a considerable challenge. In this Letter, we suggest a promising approach to create a pseudomagnetic field which induces a pure valley edge current and snake orbit in a twisted graphene nanoribbon (GNR). Such a twisting 0031-9007=14=112(9)=096805(5)
technique was successfully used in carbon nanotubes before [29,30]. The pseudomagnetic field generated in a twisted GNR can be derived by the gauge representation of electronic states [17–19,31,32]. Figure 1(a) displays a twisted GNR and its horizontal projection. It has a nearly quadratic spatial distribution of the normal strain ϵyy along the width dimension, x, Fig. 1(b). Therefore, the corresponding pseudomagnetic field for electronic states at the K valley of a GNR with chirality χ at twist rate γ is (see the Supplemental Material [33] for detailed information), BK ¼ C0 ððcos2 χ − 3sin2 χÞ cos χÞγ 2 x.
(1)
The coefficient C0 ¼ cβt0 ϕ0 a−1 0 , where t0 and a0 are the nearest neighbor hopping and the equilibrium C-C bond length of graphene, respectively. c and β [17] are both dimensionless parameters. ϕ0 is the magnetic flux quantum. It is easy to see that BK is anisotropic [34], being strongest at χ ¼ 0° (zigzag GNR), gradually weakening at 0° < χ < 30° (chiral GNR), and completely vanishing at χ ¼ 30° (armchair GNR), Fig. 1(c). On the other side, BK also varies along the width, x and changes sign at the GNR center (x ¼ 0). The pseudomagnetic field picture is instructive in elucidating the physical origin of the quantum Hall effect (QHE) and other electromagnetic properties of deformed graphene [35]. Notice that Eq. (1) is based on the nearest neighbor approximation π-orbital tight-binding model for weakly deformed graphene systems. However, generating a strong homogeneous pseudomagnetic field over a wide area usually involves a complex strain pattern and high strain level. Consequently, the dependence of π-electron hoppings between the nearest neighbors on the strain tensor ϵij ði; j ¼ x; yÞ is no longer linear, where the coupling of the electronic structure with
096805-1
© 2014 American Physical Society
PRL 112, 096805 (2014)
PHYSICAL REVIEW LETTERS
FIG. 1 (color online). (a) A twisted GNR (upper panel) and its horizontal projection (lower panel). The coordinate system is defined as x along the GNR width dimension and y along the GNR axis. Twisting is around the GNR axis. The twist-induced pseudomagnetic field of K valley-specified electronic states, BK directs reversely on the two sides of GNR axis. The pseudomagnetic field of K 0 valley electronic states BK0 ¼ −BK by time reversal symmetry. The edge currents JK and JK0 arising from K valley and K 0 valley, respectively, are counterpropagating. The current SK and SK0 associated with the snake orbit are also counterpropagating but in a serpentine manner around the GNR axis. (b) Strain distribution ϵyy along the GNR width, x at twist rate γ ¼ 0.61°=nm and (c) the corresponding pseudomagnetic field BK for 170 nm wide GNRs with different chiralities χ according to Eq. (1) with a0 ¼ 1.42 Å, β ¼ 2, c ¼ 1, and t0 ¼ 2.7 eV. (d) A near-field scanning optical microscope setup including contacts I (injector) and C (collector) for the detection of localized edge states (solid curve) and snake states (dashed curve).
elastic strain is not well defined [36,37]. Furthermore, not only nearest neighbor interactions determine the electronic structure, but also far-reaching interactions play a nonnegligible role [38]. These issues call for the investigation of the electronic states beyond the scope of the low-energy k · p Hamiltonian [2]. Atomistic quantum mechanical simulations represent a reliable alternative. Because the graphene systems used in QHE study are often of the dimension greater than 102 nm [3,17,21], having a few million atoms, full first principle and systematic quantum mechanical (QM) calculations formulated with periodic boundary conditions are impossible given the current computational capacity. Here, we showcase a study of strain-induced Landau quantization in graphene by realistic QM atomistic simulation of twisted GNRs with a width greater than 102 nm. The simulation is carried out with a density-functional theory based tight-binding (DFTB) theory [39,40]. Such a QM study is possible due to a symmetry-adapted scheme for the electronic systems [41]. Through this scheme, one can utilize the helical symmetry created by the twist to enable the realistic QM modeling of twisted GNR on a
week ending 7 MARCH 2014
basic repeating motif which contains the same N atoms in the primitive cell of the flat GNR. For example, for a zigzag GNR with a width of 170 nm, N ¼ 1608 (see Supplemental Material [33]). We have carried out DFTB calculations on a series of twisted GNRs with two types (zigzag and armchair) of edges. We optimize GNR atomistic geometries via a conjugate gradient energy minimization to the repeating motif. In our calculation, the greatest width of the GNR is about 176 nm comparable with the size of the sample used in QHE experiments [4,21]. Our result reveals a quite different picture of the Landau quantization in twisted GNRs from the one predicted by theory [17–19,31,32] as shown in Eq. (1). That is, the twist induces a well-defined electronic Landau spectrum not only in zigzag GNR but also in armchair GNR and, thus, forms localized edge states. Meanwhile, a snake orbit presents at the GNR center due to the opposite directing pseudomagnetic field on the two sides of GNR. Therefore, if any bias electrical potential is applied along the GNR axis, one shall be able to observe valley-specified edge currents, JK (JK0 ) and snake orbit currents, SK (SK0 ), schematically shown in Fig. 1(a). Figure 2 displays the band structure and density of states (DOS) of a twisted zigzag GNR of width 170 nm at (a) γ ¼ 0°=nm, and (b) γ ¼ 0.61°=nm. A very wide GNR mimics closely the band structure of bulk graphene, showing neatly linear dispersion for the valley-specified electronic states, Fig. 2(a). The extra bands associated with the zigzag edges, being dispersion free in the simple tight binding modeling with the nearest neighbor approximation, are now dispersive without the particle-hole symmetry. Mostly, this owes to the interaction between the next nearest neighbor carbon atoms [38]. Consequently, their contribution to the DOS is off the Dirac point, as indicated by the down arrow. Under twisting, we witness the onset of the flattening of valley-specified electronic bands and the DOS curve starts to fluctuate (see the online supplementary material). Thoroughly flatting of valley electronic bands starts at γ ¼ 0.61°=nm, producing sharp and separated peaks in the DOS curve, Fig. 2(b). These peaks correspond to the discrete Landau levels (LLs). It should be noted that the zero-energy LL at the Dirac point is formed from valley-specified electronic states, not from the electronic states arising from zigzag edges which, located between the zero and the first negative LLs in the DOS curve as indicated by down arrow, are not sensitive to twist. This is because they are mostly confined at very narrow regions around the GNR margins (see the Supplemental Material [33]). The twist-induced pseudomagnetic field only perturbs them weakly. Further twisting the GNR would not bring significant change in the DOS (see the Supplemental Material [33]). This means that there is a critical twist γ c , or equivalently, a saturated effective pseudomagnetic field, Beff . For this exemplified GNR, γ c ¼ 0.61°=nm. Thus, for γ > γ c , the obtained Landau spectrum stays stationary. More insight can be obtained from the fact that the effect of pseudomagnetic field on electronic states is in similar
096805-2
PRL 112, 096805 (2014)
PHYSICAL REVIEW LETTERS
FIG. 2 (color online). Electronic band structures (upper panel) and density of states (lower panel) of a 170 nm-wide zigzag GNR at twist rate (a) γ ¼ 0°=nm and (b) γ ¼ 0.61°=nm. The down arrows indicate the extra electronic states arising from zigzag edges. In (b), the solid dot (upper panel) indicates the snake orbit states, and the vertical dashed lines (lower panel) indicate the Landau levels with n ¼ 0, 1, 2, 3, 4 according to Eq. (2) with Beff ¼ 85 T. (c) The first Landau level energy, E1 , and (d) the effective magnetic field, Beff , versus GNR width W 0 .
week ending 7 MARCH 2014
A careful analysis of the experimental conditions reveals that the graphene bubbling is of a relatively small size of about only 10 nm. Over such a small size scale, it is possible to achieve a large strain gradient without requiring extremely high strain. This is consistent with our observation of larger Beff in narrower GNRs. Next we turn to the study of an armchair GNR. The theoretical prediction [17–19,31,32], as shown in Eq. (1) for twisted armchair GNRs, suggests that there would not be any Landau quantization since BK ¼ 0. Here, our result represents an exception of theory. For the flat case (γ ¼ 0°=nm), an armchair GNR displays the typical linear dispersion of electronic band structure and smooth DOS, Fig. 3(a). Without the presence of a zigzag edge, there are not any extra bands. Under twisting, the electronic states undergo the Landau quantization. Figure 3(b) shows the band structure and DOS at twist rate γ c ¼ 0.66°=nm. Again, the distribution of major DOS peaks agrees well with Eq. (2) with Beff ¼ 72 T. The nonvanishing pseudomagnetic field mainly comes from the distortion of the hexagonal lattice. Note that by time reversal symmetry, the twist-induced Landau states at K 0 valley is degenerate with those at K valley. Therefore, the corresponding pseudomagnetic field is simply −Beff . The large split of zero LL is related to the chirality property. There is also a small splitting of zero LL in twisted zigzag GNRs due to symmetry breaking of the A-B lattice. See the Supplemental Material [33] for an analysis of the symmetry breaking of the A-B lattice. A GNR mimics a Hall bar [3,21]. Because of the finite size in the width dimension, the valley-specified electronic states are pushed against the GNR margins given the strong pseudomagnetic field, Beff , forming the so-called edge
with that of a real magnetic field applied perpendicular to the graphene sheet. The formed Landau spectrum is [3,21] En ¼ sgnðnÞvF
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2eℏBeff jnj;
(2)
where the integer n indexes LLs. Indeed, with Beff ¼ 85 T, the major DOS peaks fall onto the vertical dashed lines representing the first few LLs defined by Eq. (2), Fig. 2(b). For a systematic study, we have carried out DFTB calculations on a series of GNRs with different widths. For each GNR, we determine γ c and Beff according to Eq. (2). For the evolution of the DOS with width, see the Supplemental Material [33]. We characterize the size dependence of the lowest LL and Beff by a simple fitting of the atomistic data, as E1 ∝ 1=W 0 , Fig. 2(c), and Beff ∝ 1=W 20 , Fig. 2(d), respectively. In the previous theoretical studies [17–19], the strain induced pseudomagnetic field is of tens of tesla. We suggest that this is mainly because the strain level is relatively low in these theoretical studies with the maximum strain ∼10%. Here, Beff is much larger, especially for narrow GNRs, Fig. 2(d), close to the experimental observation [20] where the strain (bubbling) induced pseudomagnetic field can be as large as 300 T.
FIG. 3 (color online). Electronic band structures (upper panel) and density of states (lower panel) of a 176 nm-wide armchair GNR at twist rate (a) γ ¼ 0°=nm and (b) γ ¼ 0.66°=nm. In (b), the solid dot (upper panel) indicates the snake states, and the vertical dashed lines (lower panel) indicate the LLs n ¼ 1, 2, 3, 4 according to Eq. (2) with Beff ¼ 136 T. The down arrows indicate the zero LL.
096805-3
PRL 112, 096805 (2014)
PHYSICAL REVIEW LETTERS
week ending 7 MARCH 2014
FIG. 5 (color online). Wave function spatial distributions (jΨj2 ) along armchair GNR width x (a) for the first four LLs n ¼ 0ðzero LLÞ, 1, 2, 3 as shown in Fig. 3(b) at κ ¼ π and (b) for the snake orbit whose location in the energy spectrum is indicated in Fig. 3(b) at κ ¼ 0.2π. In (a), the electronic states represented by the red dashed line and blue solid line are degenerate. The vertical dashed lines indicate the maxima of the wave functions for the different LLs.
FIG. 4 (color online). Wave function spatial distributions (jΨj2 ) along the zigzag GNR width x (a) for the first four LLs n ¼ 0, 1, 2, 3 as shown in Fig. 2(b) at κ ¼ 0 and (b) for the snake orbit whose location in the energy spectrum is indicated in Fig. 2(b) at κ ¼ 0.7π. In (a), the electronic states represented by the red dashed line and blue solid line are degenerate. The vertical dashed lines indicate the maxima of the wave functions for the different LLs.
states [42]. Figure 4(a) and Fig. 5(a) display the spatial distribution of wave functions of edge states along the GRN width for the studied zigzag GNR and armchair GNR, respectively. Indeed, these edge states are localized in the narrow area at the GNR margins and have different distributions for different LLs. Thus, for lower order (higher order) LLs, the valley edge currents transport via narrower (wider) atomic channels. The presence of a zigzag edge as in zigzag and chiral GNRs also brings extra edge states (see Supplemental Material [33]) which, however, are susceptible to edge disorder such as defect and reconstruction. The effects, such as elastic backscattering, are the consequences. The edge states completely insusceptible of edge disorder can be obtained in twisted armchair GNRs. This is because these states absolutely originate from valley-specified electronic states since there is no longer a zigzag edge presented, thus forming pure valley edge currents. On the other side, the oppositely directing pseudomagnetic field, Beff , on the two sides of the GNR results in a peculiar state—a snake orbit [43–45]— which is otherwise nonexistent in bent GNRs [18,19]. Spatially, it is located at the GNR center having a different momentum from edge states. See Fig. 4(b) and Fig. 5(b) for
the snake orbits in zigzag and armchair GNRs, respectively. To visualize such edge states and the snake state, a detector can be proposed utilizing state-of-art near-field scanning microwave microscopy [46,47], as schematically shown in Fig. 1(d). The twist strain can be applied in a similar way of the carbon nanotube pedal experiments [29,30]. The GNR is placed on a substrate of the insulating polymer. At the emergence of the edge current (snake current), which can be controlled by tuning the momentum of valley polarized electron beam, the highlight of the microwave image will occur at the margins (at the center) of the GNR. To summarize, we accomplished the realistic QM atomistic study of strain-induced electronic Landau quantization in GNRs under twisting, a fundamental type of deformation not considered before. Although the low energy theory [17–19,31,32] predicts that twisting is not an ideal candidate to induce Landau quantization since the strain induced pseudomagnetic field is essentially inhomogeneous (In particular it is zero for armchair GNR) as shown in Fig. 1(c), our calculation reveals that a well-defined Landau spectrum can be generated for both types (zigzag and armchair) of GNRs by twisting, demonstrating the importance of realistic atomistic calculation. We then immediately conclude that Landau quantization can also be achieved in other chiral GNRs. This revelation alleviates greatly the concern in experiment of controlling GNR’s orientation (chirality). The effective pseudomagnetic field orienting oppositely on the two sides of GNR, approaches hundreds of tesla. The electronic states, localized either at the GNRs edges, undergo Landau quantization, forming a valley-specified edge current, or at the GNR center, forming a valleyspecified snake orbit current. Note that in armchair GNR, we have pure valley edge currents. A doping effect [44] can
096805-4
PRL 112, 096805 (2014)
PHYSICAL REVIEW LETTERS
also be induced by twisting. Altogether, our results demonstrate a viable manner to control the electronic levels and transport properties of graphene at nanoscale, shedding new light on the design of novel electronics devices by strain engineering in graphene and other two-dimensional atomic crystals, e.g., silicene, MoS2 , etc. D.-B. Z. thanks Traian Dumitrica for discussion at the early stage of this work. G. S. acknowledges the support of the ERC project INTIF 226639. K. C. is supported by NSFC No. 10934007 and the National Basic Research Program of China (973 Program) under Grant No. 2011CB922204. Synergetic Innovation Center of Quantum Information and Quantum Physics at USTC, and CAS Center for excellence in quantum information and quantum physics.
*
[1]
[2] [3]
[4]
[5] [6] [7] [8] [9]
[10] [11] [12] [13] [14] [15] [16] [17] [18] [19]
Corresponding author. [email protected] 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). A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, Rev. Mod. Phys. 81, 109 (2009). K. S. Novoselov, Z. Jiang, Y. Zhang, S. V. Morozov, H. L. Stormer, U. Zeitler, J. C. Maan, G. S. Boebinger, P. Kim, and A. K. Geim, Science 315, 1379 (2007). K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A. Firsov, Science 306, 666 (2004). D. Xiao, W. Yao, and Q. Niu, Phys. Rev. Lett. 99, 236809 (2007). A. Rycerz, J. Tworzydo, and C. W. J. Beenakker, Nat. Phys. 3, 172 (2007). V. M. Pereira and A. H. Castro Neto, Phys. Rev. Lett. 103, 046801 (2009). M. A. H. Vozmediano, M. I. Katsnelson, and F. Guinea, Phys. Rep. 496, 109 (2010). S. V. Morozov, K. S. Novoselov, M. I. Katsnelson, F. Schedin, L. A. Ponomarenko, D. Jiang, and A. K. Geim, Phys. Rev. Lett. 97, 016801 (2006). M. M. Fogler, F. Guinea, and M. I. Katsnelson, Phys. Rev. Lett. 101, 226804 (2008). K. Sasaki, Y. Kawazoe, and R. Saito, Prog. Theor. Phys. 113, 463 (2005). A. F. Morpurgo and F. Guinea, Phys. Rev. Lett. 97, 196804 (2006). F. de Juan, A. Cortijo, and M. A. H. Vozmediano, Phys. Rev. B 76, 165409 (2007). F. Guinea, M. I. Katsnelson, and M. A. H. Vozmediano, Phys. Rev. B 77, 075422 (2008). E.-A. Kim and A. H. Castro Neto, Europhys. Lett. 84, 57007 (2008). F. Guinea, B. Horovitz, and P. Le Doussal, Phys. Rev. B 77, 205421 (2008). F. Guinea, M. I. Katsnelson, and A. K. Geim, Nat. Phys. 6, 30 (2009). F. Guinea, A. K. Geim, M. I. Katsnelson, and K. S. Novoselov, Phys. Rev. B 81, 035408 (2010). T. Low and F. Guinea, Nano Lett. 10, 3551 (2010).
week ending 7 MARCH 2014
[20] N. Levy, S. A. Burke, K. L. Meaker, M. Panlasigui, A. Zettl, F. Guinea, A. H. Castro Neto, and M. F. Crommie, Science 329, 544 (2010). [21] Y. Zhang, Y. -W. Tan, H. L. Stormer, and P. Kim, Nature (London) 438, 201 (2005). [22] D. A. Abanin, K. S. Novoselov, U. Zeitler, P. A. Lee, A. K. Geim, and L. S. Levitov, Phys. Rev. Lett. 98, 196806 (2007). [23] K. I. Bolotin, F. Ghahari, M. D. Shulman, H. L. Stormer, and P. Kim, Nature (London) 462, 196 (2009). [24] M. C. Rechtsman, J. M. Zeuner, A. Tünnermann, S. Nolte, M. Segev, and A. Szameit, Nat. Photonics 7, 153 (2012). [25] C. L. Kane and E. J. Mele, Phys. Rev. Lett. 95, 226801 (2005). [26] C. Lee, X. Wei, J. W. Kysar, and J. Hone, Science 321, 385 (2008). [27] K. S. Kim, Y. Zhao, H. Jang, S. Y. Lee, J. M. Kim, K. S. Kim, J. -H. Ahn, P. Kim, J. -Y. Choi, and B. H. Hong, Nature (London) 457, 706 (2009). [28] Z. H. Wu, F. Zhai, F. M. Peeters, H. Q. Xu, and K. Chang, Phys. Rev. Lett. 106, 176802 (2011). [29] T. Cohen-Karni, L. Segev, O. Srur-Lavi, S. R. Cohen, and E. Joselevich, Nat. Nanotechnol. 1, 36 (2006). [30] A. M. Fennimore, T. D. Yuzvinsky, W. -Q. Han, M. S. Fuhrer, J. Cumings, and A. Zettl, Nature (London) 424, 408 (2003). [31] J. L. Mañes, Phys. Rev. B 76, 045430 (2007). [32] H. Suzuura and T. Ando, Phys. Rev. B 65, 235412 (2002). [33] See Supplemental Material at http://link.aps.org/ supplemental/10.1103/PhysRevLett.112.096805 for the calculation method of twisted GNRs and detailed information about the spatial distributions of edge and snake orbit states. [34] T. O. Wehling, A. V. Balatsky, A. M. Tsvelik, M. I. Katsnelson, and A. I. Lichtenstein, Europhys. Lett. 84, 17003 (2008). [35] F. de Juan, A. Cortijo, M. A. H. Vozmediano, and A. Cano, Nat. Phys. 7, 810 (2011). [36] A. L. Kitt, V. M. Pereira, A. K. Swan, and B. B. Goldberg, Phys. Rev. B 85, 115432 (2012). [37] Y. Chang, T. Albash, and S. Haas, Phys. Rev. B 86, 125402 (2012). [38] N. M. R. Peres, F. Guinea, and A. H. Castro Neto, Phys. Rev. B 73, 125411 (2006). [39] D. Porezag, Th. Frauenheim, Th. Kohler, G. Seifert, and R. Kaschner, Phys. Rev. B 51, 12 947 (1995). [40] R. Ruralia and E. Hernández, Comput. Mater. Sci. 28, 85 (2003). [41] This special method was pioneered by several authors: C. T. White, D. H. Robertson, and J. W. Mintmire, Phys. Rev. B 47, 5485 (1993); V. N. Popov, New J. Phys. 6, 17 (2004); P. B. Allen, Nano Lett. 7, 1220 (2007). Also see Supplemental Material [33] for more details. [42] D. A. Abanin, K. S. Novoselov, U. Zeitler, P. A. Lee, A. K. Geim, and L. S. Levitov, Phys. Rev. Lett. 98, 196806 (2007). [43] J. E. Müller, Phys. Rev. Lett. 68, 385 (1992). [44] T. Low, F. Guinea, and M. I. Katsnelson, Phys. Rev. B 83, 195436 (2011). [45] J. R. Williams and C. M. Marcus, Phys. Rev. Lett. 107, 046602 (2011). [46] K. Lai, M. Nakamura, W. Kundhikanjana, M. Kawasaki, Y. Tokura, M. A. Kelly, and Z. -X. Shen, Science 329, 190 (2010). [47] K. Lai, W. Kundhikanjana, M. A. Kelly, Z. -X. Shen, J. Shabani, and M. Shayegan, Phys. Rev. Lett. 107, 176809 (2011).
096805-5
