Edge effects in Bilayer Graphene Nanoribbons
aa r X i v : . [ c ond - m a t . m t r l - s c i ] O c t Edge effects in Bilayer Graphene Nanoribbons
Matheus P. Lima, ∗ A. Fazzio,
1, 2, † and Antˆonio J. R. da Silva ‡ Instituto de F´ısica, Universidade de S˜ao Paulo, CP 66318, 05315-970, S˜ao Paulo, SP, Brazil. Centro de Ciˆencias Naturais e Humanas, Universidade Federal do ABC, Santo Andr´e, SP, Brazil (Dated: May 29, 2018)We show that the ground state of zigzag bilayer graphene nanoribbons is non-magnetic. It alsopossesses a finite gap, which has a non-monotonic dependence with the width as a consequence ofthe competition between bulk and strongly attractive edge interactions. All results were obtainedusing ab initio total energy density functional theory calculations with the inclusion of parametrizedvan der Waals interactions.
PACS numbers: 73.22.-f, 72.80.Rj, 61.48.De, 71.15.Nc
Since the synthesis of graphene[1], a plethora of in-triguing properties has been found in this two dimen-sional zero gap crystal due to the presence of masslessfermions with a high mobility[2, 3]. Besides the mono-layer, stacking two layers of graphene still preserves thehigh mobility, and some features of the electronic spec-trum can be controlled, for example, by applying an ex-ternal electric field[4]. Measurements of quantum halleffect and quasi-particle band structure indicate qualita-tive differences between monolayers and bi-layers. Theoccurrence of this rich physics at room-temperature[5, 6]has attracted a great interest in designing graphene-based nanoelectronic devices. In this scenario, it is fun-damental to establish and control an energy gap ( E g ).A possibility is to introduce lateral quantum confine-ment via synthesis of single layer (GNR)[7, 8] or bilayer(B-GNR)[9, 10] graphene nanoribbons by plasma etch-ing or chemical routes[11, 12]. This opens a gap thatrises the possibilities to use graphene in nanoelectronics,where small widths (sub-10 nm ) is required for room tem-perature applications[6]. However, the BGRNs are lesssensitive to external perturbations in comparison withGNRs, and hence, they may be more appropriate to fab-ricate hight quality nano-devices[13].The electronic structure of the nanoribbons, includingthe gap, are largely affected by the geometric pattern(zigzag or armchair) at their edges. GNRs with zigzagedges, in particular, have as a distinct feature the pres-ence of edge states that introduce a large density of states(DOS) at the Fermi energy. Theoretical works predictthat this configuration is unstable, and there will be theappearance of a magnetic order that leads to a removalof this large DOS peak[7, 8, 14]. Magnetism in GNRs hasbeen intensively investigated as a possible way to developspintronic devices[15, 16, 17]. An anti-ferromagnetic (fer-romagnetic) order between the two edges leads to a semi-conductor (metallic) state[7, 8, 14]. It is also believedthat magnetism is necessary to open a gap in zigzag B-GNR (B-ZGNR)[18].In this Letter we investigate the geometrical and elec-tronic structure of B-ZGNR, and show that the groundstate of these systems is non-magnetic and possesses a non-monotonic finite gap . There are two possible edgealignments for the B-ZGNR, called α and β (see Fig.1). We found that the α alignment is energetically fa-vorable, with an inter-layer edges attraction, whereasfor the β alignment there is an inter-layer edges repul-sion. These edge-related forces cause a deviation fromthe exact Bernal stacking, resulting in a non-monotonicbehavior of the energy gap with the width w for the α B-ZGNR, with a maximum value at w ≈ . nm . Theseresults differ qualitatively from GNRs with zigzag edges(M-ZGNR)[11, 12]All our results are based on ab initio total energy Den-sity Functional Theory[19] (DFT) calculations. In orderto correctly describe multi-layer graphitic compounds, itis necessary to include van der Waals (vdW) interactions.The use of fully relaxed total energy DFT calculations tostudy such systems suffers from serious limitations, sincethe most traditional exchange correlation ( xc ) function-als in use today do not correctly describe these terms.With the LDA xc , the geometry is correctly describedbut the inter-layer binding energy is underestimated by hd d ( ) a ( ) b ( ) c FIG. 1: Bilayer graphene nanoribbons with (a) α and (b) β edge alignment. (c) Side view of bilayer graphene nanorib-bons. The blue (white) atoms form the upper (bottom) layer. TABLE I: Inter-layer binding energy ( E b ), in eV /atom , anddistance h , in ˚A, for graphite and a graphene bilayer, whichis not bound (NB) at the GGA level.present Exp LDA GGA(PBE)Graphite E b h 0.0543.350 0 . ± . a . b E b h 0.0273.320 -- 0.0173.202 NBNB a From Ref. [21] b From Ref. [22] xc do not even correctly de-scribe the geometrical features[20]. Thus, in order to beable to investigate the geometries and relative energies ofB-ZGNRs, we include a non-local potential in the Kohn-Sham (KS) equations that correctly describes the vdWinteractions. We modified the SIESTA code[23] addingin the KS Hamiltonian[24] the dispersion corrected atomcentered potential (DCACP)[25]. This correction is suf-ficiently accurate to describe weakly bonded systems[26]with the vdW interactions included in the whole self-consist cycle, providing accurate values for both forces(and thus geometries) and total energies. Our imple-mentation was successfully tested (see Table I), and wasemployed to obtain the results here reported[28, 29, 30].We investigated B-ZGNR composed by two M-ZGNRpassivated with hydrogen, and widths[31] that rangefrom w = 0 . w = 4 . nm . The layers are in theBernal stacking, which means that there are two types ofC atoms, those that are positioned above the center ofthe hexagons of the other layer, defining a B-sublattice,and those right on top of the C atoms of the other layer,forming an A-sublattice. An infinite graphene bilayer hasno gap, and the orbitals at the Fermi level are located atthe B-sublattice. When we cut the layer along the zigzagedge, there are two possible alignments (Fig.1): (a) the α alignment, where the outermost edge atoms belong to theA-sublattice, and (b) the β alignment, where the outer-most edge atoms belong to the B-sublattice. Thus, onlythe inter-layer edge interaction differs. Two geometri-cal distortions have proven to be important: (i) an edgedistortion that causes a curvature in the ribbons (seeFig.1(c)), and (ii) a lateral deviation from the perfectBernal stacking. To quantify this deviation, we definethe quantity u ≡ d C − C − d , where d is shown in Fig.1(a)and (b), and d C − C is the carbon-carbon bond length.The perfect Bernal stacking corresponds to u = 0.The geometries and band structures of fully relaxedB-ZGNRs with w ≈ . nm for α and β alignments arepresented in Fig. 2. In B-ZGNR with the β alignment,similarly to M-ZGNRs, a non spin-polarized calculationsleads to a high DOS at the Fermi energy, and a magnetic FIG. 2: Ground State of fully relaxed B-ZGNRs generated bystacking two (5,0) M-ZGNR. Below each band structure thegeometry and local magnetization are presented. (a) α align-ment. This state is non-magnetic and presents a geometricdistortion near the edge. (b) β alignment. This state showsan AF in-layer and AF inter-layer magnetic order. order is required to split these localized edge states atthe K symmetry point. In order to establish the pos-sible spin polarized configurations, we used four initialguesses for the density matrix before starting the self-consistency cycle, which are: i ) anti-ferromagnetic (AF)in-layer and inter-layer, ii ) ferromagnetic (F) in-layer andinter-layer, iii ) AF in-layer and F inter-layer and iv ) F in-layer and AF inter-layer. As well as non-polarized calcu-lations. From all calculations, the AF in-layer and inter-layer guess leads to the lower energy state (Fig. 2(b)).However, the energy differences are less than k B T [18].At the α alignment (Fig. 2(a)), on the other hand, weobtain a qualitatively new situation. There is a strongattractive interaction between the edge atoms of the twolayers, with a resulting geometric distortion that de-creases the distance between them (Fig. 2(a)). For all α B-ZGNRs that we have investigated, the final geometryalways had an inter-layer edge atoms distance around 3 . A . The final configuration is non-magnetic and with a fi-nite gap, contrary to previous results where the presenceof a gap was intrinsically coupled to a magnetic state[18].Note that if we do not allow the atoms at the two layers torelax, but simply optimize the inter-layer distance ( i.e. ,the layers keep their planar geometries), a magnetic con-figuration is still necessary to open a gap[18]. However,this configuration has higher energy.If we take one of the mono-layers that form the finalrelaxed α B-ZGNR, and perform a calculation withoutletting the atoms relax, we obtain an energy increase,when compared to the lowest energy M-ZGNR (AF in-layer[15]), that can be broken down in two components.Considering M-ZGNR with widths larger than 2 nm , ifwe allow the distorted monolayer to be magnetic, the en-ergy increase is ≈ . eV /nm , which can be viewed as theelastic contribution. If we now consider a non-magneticconfiguration for this distorted monolayer, which is thesituation in the α B-ZGNR, there is an extra energy in-crease of ≈ . eV /nm , i.e. , an overall energy penalty of ≈ . eV /nm . Considering the two monolayers, the totalenergy cost to deform and demagnetize the α B-ZGNR is ≈ . eV /nm . This energy increase is more than com-pensated by the edge atoms interaction, and the energygain is in part associated with a large split of the localizededge states. At the β B-ZGNR (Fig. 2(b)), the uniqueway to diminish the DOS associated with the localizededge states at the Fermi energy is via a magnetic order-ing, and the system tends to increase the inter-layer edgeatoms distance in order to allow a bigger magnetization,given rise to a repulsive inter-layer edge interaction.
FIG. 3: Dependence of the binding energies with the width.
Comparing the two alignments, the α B-ZGNR resultsenergetically favorable. This is an even more importantconclusion considering that most of the calculations usedthe β B-ZGNR[10, 32]. Fig. 3 presents the dependence ofthe B-ZGNRs binding energies, for both edge alignments,as a function of w (calculated relative to two isolated low-est energy M-ZGNR). The interaction between the layerscan be separated into two components: i ) edge interac-tions, that do not depends on the width, ii ) and bulkinteractions, that increase linearly with w . The bindingenergies (per unit length) can be well adjusted with: E b ( w ) = a + bw. (1)Since there are two edges, a/ b is the inter-layerbulk interaction energy per unit area. At the α align-ment, a = − . eV /nm and b = − . eV /nm , indicating FIG. 4: Dependence of the (a) energy gap (inset illustrates thechange of character of the VBM), and (b) the lateral deviation u with the width w (inset indicates how u affects the inter-layer A-sublattice p -orbitals interaction). that the edges interaction is attractive ( a < β alignment, a = +0 . eV /nm and b = − . eV /nm , indi-cating that there is a repulsive edges interaction ( a > b , as expected, does not depend on theedge alignment, and it is very close to our calculatedbulk inter-layer interaction in a graphene bilayer (0.027 eV /atom = − . eV /nm ).For the α B-ZGNR, there is a competition between theforces deriving from the bulk, that do prefer the Bernalpattern of stacking, and the forces deriving from theedges, that tend to maintain the inter-layer edge carbonatoms distance close to 3 ˚A. The system then minimizesthe overall energy penalty by simultaneously optimizingboth the deviation u from the exact Bernal stacking, andthe elastic energy associated with the ribbon’s curvature.Thus, for narrower B-ZGNR the system prefers to have alarger value of u and a smaller overall curvature whereasfor wider B-ZGNR the deviation u tends to decrease tominimize the bulk energy penalty, since now it is possibleto have a softer curvature that is somewhat localized atthe edges when compared to the total width of the rib-bon (see Fig. 1(c)). As a result, the dependence of thelateral deviation u with the width w is well adjusted by u ( w ) = 0 . e − . w (with u in [˚ A ], and w in [ nm ]). More-over, as a result, the carbon-carbon bond lengths do notsignificantly differ from their values in the M-ZGNRs.For this situation, the average inter-layer distance h isclose to its value in the graphene bilayer (see Tab. I).Without the geometrical deformation caused by theinter-layer edge interactions, a monotonic decrease of theenergy gap is expected due to the quantum confinement( ∝ /w )[8, 18]. However, for small ribbons, we find thatthe character of the valence band maximum (VBM) islocated at the A-sublattice, as opposed to larger rib-bons where it is located in the B-sublattice (see Fig.4(a)), similarly to the Fermi level orbitals in the infi-nite graphene bilayer. Moreover, since the interactionbetween the C atoms in the A-sublattice increases when u decreases (see Fig.4(b)), the gap initially increases with w . However, for w > ∼ . nm , due to the quantum-confinement decrease, there is a crossover between thetwo highest occupied bands, and the character of theVBM is at the B-sublattice. Thus, this leads to a non-monotonic behavior of the energy gap with w (Fig. 4(a)).For the β B-ZGNR, there is a repulsive interactionbetween the edges, in such way that the inter-layer edgecarbon atoms distance is close to 3.7 ˚A. There happensa small negative lateral deviation ( u <
0) that can beneglected when w > . nm . Note that for the α B-ZGNR, due to the attractive edges interaction, u > w > nm .Summarizing, we unequivocally show that for B-ZGNRthe edge alignment α is the lowest energy configuration.This is a result of the strong attractive interaction be-tween the edges, that is manifested in an observed chem-ical bonding between the inter-layer edge carbon atoms,and that significantly influences the geometry and elec-tronic structure of bilayer nanoribbons with sub-10 nm widths. As a consequence, the ground-state is non-magnetic and possesses a finite gap, which presents anon-monotonic dependence with the width.We acknowledge helpful discussions with M. D.Coutinho-Neto regarding the DCACP, financial supportfrom FAPESP and CNPq. ∗ [email protected] † [email protected] ‡ [email protected][1] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang,Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A.Firsov, Science , 366 (2004).[2] 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, Nature , 197 (2005); Y. B.Zhang, Y. W. Tan, H. L. Stormer, P. Kim, Nature ,201 (2005).[3] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S.Novoselov, and A. K. Geim, Rev. Mod. Phys., accepted.[4] K. S. Novoselov, E. McCann, S. V. Morozov, V. I. Fal’ko,M. I. Katsnelson, U. Zeitler, D. Jiang, F. Schedin, A. K.Geim, Nat. Phys. , 177 (2006); T. Ohta, A. Bostwick, T.Seyller, K. Horn, E. Rotenberg, Science , 951 (2006);J. B. Oostinga, H. B. Heersche, X. L. Liu, A. F. Mor-purgo, L. M. K. Vandersypen, Nat. Mat. , 151 (2008).[5] K. S. Novoselov, Z. Jiang, Y. Zhang, S. V. Morozov, H.L. Stormer, U. Zeitler, J. C. Maan, G. S. Boebinger, PKim, A. K. Geim, Science , 1379 (2007); N. Tombros,C. Jozsa, M. Popinciuc, H. T. Jonkman, B. J. van Wees,Nature , 571 (2007).[6] X. R. Wang, Y. J. Ouyang, X. L. Li, H. L. Wang, J. Guo,H. J. Dai, Phys. Rev. Lett. , 206803 (2008).[7] Y. W. Son, M. L. Cohen, S. G. Louie, Phys. Rev. Lett. , 216803 (2006). [8] L. Yang, C. H. Park, Y. W. Son, M. L. Cohen, S. G.Louie, Phys. Rev. Lett. , 186801 (2007).[9] K.-T. Lam and G. Liang, Appl. Phys. Lett. , 223106(2008).[10] E. V. Castro, N. M. R. Peres, and J. M. B. Lopes dosSantos, JOAM , 1716 (2008).[11] M. Y. Han, B. Ozyilmaz, Y. B. Zhang, P. Kim, Phys.Rev. Lett. , 206805 (2007).[12] X. L. Li, X. R. Wang, L. Zhang, S. W. Lee, H. J. Dai,Science , 1229 (2008).[13] Y. M. Lin, P. Avouris, Nano Lett. , 2119 (2008).[14] L. Pisani, J. A. Chan, B. Montanari, and N. M. Harrison,Phys. Rev. B , 064418 (2007).[15] T. B. Martins, R. H. Miwa, A. J. R. da Silva, A. Fazzio,Phys. Rev. Lett. , 196803 (2007).[16] Y. W. Son, M. L. Cohen, S. G. Louie, Nature , 347(2006)[17] E. J. Kan, Z. Y. Li, J. L. Yang, J. G. Hou, J. Am. Chem.Soc. , 4224 (2008)[18] B. Sahu, H. Min, A. H. MacDonald, S. K. Banerjee, Phys.Rev. B , 045404 (2008).[19] P. Hohenberg, W. Kohn ,Phys. Rev. , B864 (1964);W. Kohn and L.J. Sham, Phys. Rev. , A133 (1965).[20] M. Hasegawa, K. Nishidate, H. Iyetomi, Phys. Rev. B. , 115424 (2007) and references therein.[21] R. Zacharia, , H. Ulbricht, T. Hertel, Phys. Rev. B ,155406 (2004)[22] A. Bosak, M. Krisch, M. Mohr, J. Maultzsch, C. Thom-sen, Phys. Rev. B , 153408 (2007)[23] E. Artacho, D. Sanchez-Portal, P. Ordejon, A. Garcia, J.M. Soler, Phys. Status Solidi B , 809 (1999).[24] The added matrix elements are H DCACPµν = P j P lm = − l h ϕ µ p lmj i σ j h p lmj ϕ ν i , where { ϕ µ } are thelocalized basis. We use only l =3, and centeredat each atom j there is a normalized projector p lmj ( r ) ∝ Y lm (ˆ r ) r l exp (cid:16) − r σ (cid:17) , where Y lm are sphericalharmonics. The parameters ( σ , σ ) are from Ref. [27].[25] O. A. von Lilienfeld, I. Tavernelli, U. Rothlisberger, D.Sebastiani, Phys. Rev. Lett , 153004 (2004).[26] E. Tapavicza, I. C. Lin, O. A. von Lilienfeld, I. Taver-nelli, M. D. Coutinho-Neto, U. Rothlisberger, J. Chem.Theory Comput. , 1673 (2007); O. A. von Lilienfeld, I.Tavernelli, U. Rothlisberger, D. Sebastiani, Phys. Rev. B , 195119 (2005)[27] I. C. Lin, M. D. Coutinho-Neto, C. Felsenheimer, O. A.von Lilienfeld, I. Tavernelli, U. Rothlisberger, Phys. Rev.B , 205131 (2007).[28] We use a generalized gradient aproximation[29], norm-conserving pseudopotentials[30], a mesh cut off of 400 Ry for the grid integration, a 1 × ×
50 mesh for the firstBrillouin zone sampling, a force criterion of 10 meV / ˚ A ,and a supercell approximation with a lateral separationbetween the images of 20 ˚ A .[29] J. P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. , 3865 (1996).[30] N. Troullier and J. L. Martins, Phys. Rev. B , 1993(1991).[31] We have defined the width as the lateral distance betweenthe outermost C edge atoms of the bilayer.[32] E. V. Castro, N. M. R. Peres, J. M. B. Lopes dos Santos,A. H. Castro Neto, F. Guinea, Phys. Rev. Lett. ,026802 (2008). J. Rhim, K. Moon, J. Phys. Condens.Matter,20