Table Of Content

Boundary conditions forDiracfermions onaterminated honeycomb lattice A. R. Akhmerov and C. W. J. Beenakker Instituut-Lorentz, Universiteit Leiden, P.O.Box 9506, 2300 RALeiden, TheNetherlands (Dated:October,2007) WederivetheboundaryconditionfortheDiracequationcorrespondingtoatight-bindingmodelonatwo- dimensionalhoneycomblatticeterminatedalonganarbitarydirection.Zigzagboundaryconditionsresultgener- icallyonce theboundary isnot parallel tothebonds. Sinceahoneycomb stripwithzigzagedges isgapless, thisimpliesthatconfinementbylatticeterminationdoesnotingeneralproduceaninsulatingnanoribbon. We 8 considertheopeningofagapinagraphenenanoribbonbyastaggeredpotentialattheedgeandderivethecor- 0 respondingboundaryconditionfortheDiracequation.Weanalyzetheedgestatesinananoribbonforarbitrary 0 boundaryconditionsandidentifyaclassofpropagatingedgestatesthatcomplementtheknownlocalizededge 2 statesatazigzagboundary. n a PACSnumbers:73.21.Hb,73.22.Dj,73.22.-f,73.63.Bd J 1 1 I. INTRODUCTION dencecharacteristicofinsulatingbehaviorrequiresthespecial armchairorientation(ϕamultipleof60◦),atwhichthedecay l] Theelectronicpropertiesof graphenecanbe describedby ratef(ϕ)vanishes. al a differenceequation (representinga tight-bindingmodelon Confinement by a mass term in the Dirac equation does h a honeycomb lattice) or by a differential equation (the two- produceanexcitationgapregardlessoftheorientationofthe s- dimensionalDiracequation)[1,2]. Thetwodescriptionsare boundary. We show how the infinite-mass boundary condi- e equivalent at large length scales and low energies, provided tion of Ref. [8] can be approached starting from the zigzag m the Dirac equation is supplemented by boundary conditions boundary condition, by introducing a local potential differ- . consistentwiththetight-bindingmodel.Theseboundarycon- enceonthetwo sublatticesin thetight-bindingmodel. Such t a ditionsdependon a variety of microscopicproperties, deter- a staggered potential follows from atomistic calculations[3] m minedbyatomisticcalculations[3]. andmaywellbetheoriginoftheinsulatingbehaviorobserved experimentallyingraphenenanoribbons[9,10]. - For a general theoretical description, it is useful to know d whatboundaryconditionsonthe Dirac equationare allowed Theoutlineofthispaperisasfollows.InSec.IIweformu- n by the basic physical principles of current conservation and late,followingRefs.[4,5],thegeneralboundaryconditionof o (presence or absence of) time reversal symmetry — inde- theDiracequationonwhichouranalysisisbased. InSec.III c [ pendently of any specific microscopic input. This problem we derivefrom the tight-bindingmodelthe boundarycondi- was solved in Refs. [4, 5]. The general boundary condi- tioncorrespondingtoanarbitrarydirectionoflatticetermina- 3 tion depends on one mixing angle Λ (which vanishes if the tion.InSec.IVweanalyzetheeffectofastaggeredboundary v 3 boundarydoesnotbreaktimereversalsymmetry),onethree- potentialon the boundarycondition. In Sec. V we calculate 2 dimensionalunitvectornperpendiculartothenormaltothe the dispersion relation for a graphene nanoribbon with arbi- 7 boundary, and one three-dimensional unit vector ν on the trary boundary conditions. We identify dispersive (= propa- 2 Blochsphereofvalleyisospins. Altogether,fourrealparame- gating)edgestateswhichgeneralizetheknowndispersionless 0. tersfixtheboundarycondition. (= localized)edge states at a zigzag boundary[11]. The ex- 1 Inthepresentpaperweinvestigatehowtheboundarycondi- ponentialdependenceof the gap∆ onthe nanoribbonwidth 7 tiondependsonthecrystallographicorientationofthebound- iscalculatedinSec.VIbothanalyticallyandnumerically.We 0 ary. As the orientation is incremented by 30◦ the boundary concludeinSec.VII. : v configurationswitchesfromarmchair(paralleltoone-thirdof i thecarbon-carbonbonds)tozigzag(perpendiculartoanother X one-thirdofthebonds).Theboundaryconditionsforthearm- II. GENERALBOUNDARYCONDITION r chair and zigzag orientations are known [6]. Here we show a that the boundarycondition for intermediate orientationsre- Thelong-wavelengthandlow-energyelectronicexcitations mainsofthezigzagform,sothatthearmchairboundarycon- ingraphenearedescribedbytheDiracequation ditionisonlyreachedforadiscretesetoforientations. Sincethezigzagboundaryconditiondoesnotopenupagap HΨ=εΨ (2.1) in the excitation spectrum [6], the implication of our result (notnoticedinearlierstudies[7])isthataterminatedhoney- withHamiltonian comblatticeofarbitraryorientationismetallicratherthanin- sulating.Wepresenttight-bindingmodelcalculationstoshow H =vτ (σ p) (2.2) 0 ⊗ · that, indeed, the gap∆ exp[ f(ϕ)W/a] in a nanoribbon ∝ − atcrystallographicorientationϕvanishesexponentiallywhen actingona four-componentspinorwave functionΨ. Here v its width W becomes large compared to the lattice constant is the Fermi velocity and p = i~ is the momentum op- − ∇ a,characteristicofmetallicbehavior. The∆ 1/W depen- erator. Matricesτ ,σ are Paulimatricesin valleyspace and i i ∝ 2 sublatticespace,respectively(withunitmatricesτ ,σ ). The III. LATTICETERMINATIONBOUNDARY 0 0 currentoperatorinthedirectionnisn J =vτ (σ n). 0 · ⊗ · The Hamiltonian H is written in the valley isotropic rep- Thehoneycomblatticeofacarbonmonolayerisatriangu- resentation of Ref. [5]. The alternative representationH′ = larlattice(latticeconstanta)withtwoatomsperunitcell,re- vτz (σ p)ofRef. [4]is obtainedbytheunitarytransfor- ferredtoasAandBatoms(seeFig.1a). TheAandBatoms ⊗ · mation separately form two triangular sublattices. The A atoms are connectedonlytoBatoms,andviceversa. Thetight-binding H′ =UHU†, U = 1(τ +τ ) σ +1(τ τ ) σ . (2.3) 2 0 z ⊗ 0 2 0− z ⊗ z equationsonthehoneycomblatticearegivenby As described in Ref. [4], the general energy-independent εψ (r)=t[ψ (r)+ψ (r R )+ψ (r R )], boundaryconditionhastheformofalocallinearrestrictionon A B B 1 B 2 − − (3.1) thecomponentsofthespinorwavefunctionattheboundary: εψB(r)=t[ψA(r)+ψA(r+R1)+ψA(r+R2)]. Ψ=MΨ. (2.4) Here t is the hopping energy, ψ (r) and ψ (r) are the A B electron wave functions on A and B atoms belonging to The 4 4 matrix M has eigenvalue 1 in a two-dimensional × the same unit cell at a discrete coordinate r, while R1 = subspace containing Ψ, and without loss of generality we (a√3/2, a/2), R = (a√3/2,a/2) are lattice vectors as 2 mayassumethatM haseigenvalue 1intheorthogonaltwo- − showninFig.1a. − dimensionalsubspace. ThismeansthatM maybechosenas Regardlessofhowthelatticeisterminated,Eq.(3.1)hasthe aHermitianandunitarymatrix, electron-holesymmetryψ ψ ,ε ε. Forthelong- B B → − → − M =M†, M2 =1. (2.5) wavelength Dirac Hamiltonian (2.2) this symmetry is trans- latedintotheanticommutationrelation Therequirementofabsenceofcurrentnormaltothebound- ary, Hσz τz +σz τzH =0. (3.2) ⊗ ⊗ Ψn J Ψ =0, (2.6) Electron-holesymmetryfurtherrestrictstheboundarymatrix B h | · | i M inEq.(2.10)totwoclasses: zigzag-like(ν = zˆ,n=zˆ) with nB a unit vector normal to the boundary and pointing and armchair-like (ν = n = 0). In this sect±ion we will z z outwards,isequivalenttotherequirementofanticommutation showthatthe zigzag-likeboundaryconditionappliesgeneri- ofthematrixM withthecurrentoperator, callytoanarbitraryorientationofthelatticetermination.The armchair-likeboundaryconditionis only reachedfor special M,n J =0. (2.7) { B · } orientations. ThatEq.(2.7)impliesEq.(2.6)followsfrom Ψn J Ψ = B h | · | i ΨM(n J)M Ψ = Ψn J Ψ . The converse is B B hpro|veninA·pp.A.| i −h | · | i A. Characterizationoftheboundary Wearenowfacedwiththeproblemofdeterminingthemost general4 4matrixM thatsatisfiesEqs.(2.5)and(2.7). Ref. Aterminatedhoneycomblatticeconsistsofsiteswiththree × [4]obtainedtwofamiliesoftwo-parametersolutionsandtwo neighborsintheinteriorandsiteswithonlyoneortwoneigh- more families of three-parameter solutions. These solutions bors at the boundary. The absent neighboring sites are in- are subsets of the single four-parameter family of solutions dicated by open circles in Fig. 1 and the dangling bonds by obtainedinRef.[5], thin line segments. The tight-binding model demands that the wave function vanishes on the set of absent sites, so the M =sinΛτ (n σ)+cosΛ(ν τ) (n σ), (2.8) 0 1 2 firststepinouranalysisisthecharacterizationofthisset. We ⊗ · · ⊗ · assume that the absent sites form a one-dimensional super- whereν,n ,n arethree-dimensionalunitvectors,suchthat 1 2 lattice, consisting ofa supercellof N emptysites, translated n and n are mutually orthogonal and also orthogonal to 1 2 overmultiplesofasuperlatticevectorT. Sincetheboundary n . Aproofthat(2.8)isindeedthemostgeneralsolutionis B superlattice is part of the honeycomb lattice, we may write giveninApp.A. OnecanalsocheckthatthesolutionsofRef. [4]aresubsetsofM′ =UMU†. T = nR1 +mR2 withnandm non-negativeintegers. For example, in Fig. 1 we have n = 1, m = 4. Without loss In this work we will restrict ourselvesto boundarycondi- ofgenerality,andforlaterconvenience,wemayassumethat tionsthatdonotbreaktimereversalsymmetry. Thetimere- m n=0(modulo3). versaloperatorinthevalleyisotropicrepresentationis − The angle ϕ between T and the armchair orientation (the T = (τ σ ) , (2.9) x-axisinFig.1)isgivenby y y − ⊗ C with the operator of complex conjugation. The boundary 1 n m π π C ϕ=arctan − , ϕ . (3.3) condition preserves time reversal symmetry if M commutes (cid:18)√3n+m(cid:19) −6 ≤ ≤ 6 withT. ThisimpliesthatthemixingangleΛ = 0,sothatM isrestrictedtoathree-parameterfamily, The armchair orientation correspondsto ϕ = 0, while ϕ = π/6correspondstothezigzagorientation. (Becauseofthe M =(ν τ) (n σ), n nB. (2.10) ±π/3periodicityweonlyneedtoconsider ϕ π/6.) · ⊗ · ⊥ | |≤ 3 with ~k = p T. While the continuous quantum number · k (0,2π) describes the propagation along the boundary, ∈ a second (discrete) quantum number λ describes how these boundarymodesdecayawayfromtheboundary. Weselectλ bydemandingthattheBlochwave(3.4)isalsoasolutionof ψ(r+R )=λψ(r). (3.5) 3 ThelatticevectorR = R R hasanonzerocomponent 3 1 2 − acosϕ > a√3/2 perpendicular to T. We need λ 1 to | | ≤ preventψ(r)fromdivergingintheinteriorofthelattice. The decaylengthl inthedirectionperpendiculartoT isgiven decay by acosϕ l = − . (3.6) decay ln λ | | FIG. 1: (a) Honeycomb latice constructed from a unit cell (grey TheboundarymodessatisfyingEqs.(3.4)and(3.5)arecal- rhombus) containing twoatoms(labeledAandB), translatedover culatedin App.B fromthe tight-bindingmodel. In the low- latticevectors R1 and R2. Panels b,c,d show three different peri- energy regime of interest (energies ε small compared to t) odicboundarieswiththesameperiodT =nR1+mR2.Atomson there is an independentset of modeson each sublattice. On the boundary (connected by thick solidlines) have dangling bonds sublatticeAthequantumnumbersλandkarerelatedby (thindottedlinesegments)toemptyneighboringsites(opencircles). The number N of missing sites and N′ of dangling bonds per pe- ( 1 λ)m+n =exp(ik)λn (3.7a) riod is ≥ n+m. Panel d shows a minimal boundary, for which − − N =N′ =n+m. andonsublatticeBtheyarerelatedby ( 1 λ)m+n =exp(ik)λm. (3.7b) The number N of empty sites per period T can be arbi- − − trarily large, but it cannot be smaller than n + m. Like- For a given k there are roots λ of Eq. (3.7a) having A p wise, the number N′ of dangling bonds per period cannot absolutevalue 1,withcNorrespondingboundarymodesψp. ≤ be smaller than n + m. We call the boundary minimal if We sort these modes according to their decay lengths from N = N′ = n+m. For example, the boundary in Fig. 1d short to long, ldecay(λp) ldecay(λp+1), or λp λp+1 . ≤ | | ≤ | | ′ ThewavefunctiononsublatticeAisasuperpositionofthese isminimal(N = N = 5), whilethe boundariesin Figs.1b ′ ′ modes and1carenotminimal(N =7,N =9andN =5,N =7, respectively). In what follows we will restrict our consider- NA ations to minimal boundaries, both for reasons of analytical ψ(A) = α ψ , (3.8) p p simplicity[12]andforphysicalreasons(itisnaturaltoexpect Xp=1 thattheminimalboundaryisenergeticallymostfavorablefor agivenorientation). withcoefficientsα suchthatψ(A) vanishesontheN miss- p A ′ We conclude this subsection with a property of minimal ingAsites. Similarlythereare rootsλ ofEq.(3.7b)with boundaries that we will need later on. The N empty sites λ′ 1, λ′ λ′ . ThecNorBrespondinpgboundarymodes per period can be divided into NA empty sites on sublattice |forpm|≤thew|avpe|≤fun|ctpi+on1|onsublatticeB, AandN emptysitesonsublatticeB. Aminimalboundary B isconstructedfromntranslationsoverR ,eachcontributing NB oneemptyAsite,andmtranslationsover1R2,eachcontribut- ψ(B) = α′pψp′, (3.9) ing one empty B site. Hence, N = n and N = m for a Xp=1 A B minimalboundary. withα′ suchthatψ(B)vanishesontheN missingBsites. p B B. Boundarymodes C. Derivationoftheboundarycondition Theboundarybreaksthetwo-dimensionaltranslationalin- To derivetheboundaryconditionfortheDirac equationit varianceoverR1andR2,butaone-dimensionaltranslational is sufficient to consider the boundary modes in the k 0 invariance over T = nR1 + mR2 remains. The quasimo- limit. The characteristic equations (3.7) for k = 0→each mentum p along the boundary is therefore a good quantum have a pair of solutions λ = exp( 2iπ/3) that do not ± ± number.ThecorrespondingBlochstatesatisfies depend on n and m. Since λ = 1, these modes do ± | | not decay as one moves away from the boundary. The cor- ψ(r+T)=exp(ik)ψ(r), (3.4) responding eigenstate exp( iK r) is a plane wave with ± · 4 wave vector K = (4/3)πR /a2. One readily checks that D. Precisionoftheboundarycondition 3 this Bloch state also satisfies Eq. (3.4) with k = 0 [since K T =2π(n m)/3=0(modulo2π)]. Ataperfectzigzagorarmchairedgethefourcomponentsof · − Thewavefunctions(3.8)and(3.9)onsublatticesAandB theDiracspinorΨaresufficienttomeettheboundarycondi- inthelimitk 0taketheform tion.Neartheboundarieswithlargerperiodandmorecompli- → catedstructurethewavefunction(3.10)alsonecessarilycon- NA−2 ′ ψ(A) =Ψ eiK·r+Ψ e−iK·r+ α ψ , (3.10a) tainsseveralboundarymodesψp,ψpthatdecayawayfromthe 1 4 p p boundary.Thedecaylengthδoftheslowestdecayingmodeis Xp=1 the distance at which the boundaryis indistinguishablefrom ψ(B) =Ψ eiK·r+Ψ e−iK·r+NB−2α′ψ′. (3.10b) aperfectarmchairorzigzagedge.Atdistancessmallerthanδ 2 3 p p theboundaryconditionbreaksdown. Xp=1 Inthecaseofanarmchair-likeboundary(withn= m),all ′ thecoefficientsα andα inEqs. (3.10)mustbenonzeroto The four amplitudes(Ψ , iΨ , iΨ , Ψ ) Ψ form the p p 1 − 2 3 − 4 ≡ satisfy theboundarycondition. Themaximaldecaylengthδ four-componentspinorΨintheDiracequation(2.1). There- isthenequaltothedecaylengthoftheboundarymodeψ maining 2and 2termsdescribedecayingboundary n−1 NA− NB− whichhas the largest λ. Itcan be estimated fromthe char- modesofthetight-bindingmodelthatarenotincludedinthe | | acteristic equations(3.7) that δ T . Hence the larger the Diracequation. ≈ | | period of an armchair-like boundary, the larger the distance We are now ready to determine what restriction on Ψ is from the boundary at which the boundary condition breaks imposedby the boundaryconditionon ψ(A) and ψ(B). This down. restriction is the required boundary condition for the Dirac Forthezigzag-likeboundarythesituationis different. On equation.InApp.Bwecalculatethat,fork =0, onesublatticetherearemoreboundarymodesthanconditions imposed by the presence of the boundary and on the other =n (n m)/3+1, (3.11) NA − − sublatticetherearelessboundarymodesthanconditions. Let B =m (m n)/3+1, (3.12) usassume thatsublattice A hasmore modesthan conditions N − − (which happens if n < m). The quickest decaying set of sothat A+ B =n+m+2isthetotalnumberofunknown boundarymodessufficienttosatisfythetight-bindingbound- N N amplitudesin Eqs. (3.8) and (3.9). These have to be chosen ary condition contains n modes ψ with p n. The dis- p such that ψ(A) and ψ(B) vanish on NA and NB lattice sites tance δ from the boundary within which the≤boundary con- respectively. For the minimal boundaryunder consideration dition breaks down is then equal to the decay length of the we haveNA = n equationstodetermine A unknownsand slowestdecayingmodeψ inthissetandisgivenby N n N =mequationstodetermine unknowns. B B N Three cases can be distinguished [in each case n m = δ =l (λ )= acosϕ/ln λ . (3.13) decay n n − − | | 0(modulo3)]: [SeeEq.(3.6).] 1. If n > m then n and m+2, so Ψ = As derivedin App. B for the case of large periods T Ψ4 =0,whileΨN2Aan≤dΨ3areNunBde≥termined. 1 a, the quantum number λn satisfies the following sys|tem| ≫of equations: 2. If n < m then n and m+2, so Ψ = NB ≤ NA ≥ 2 1+λ m+n = λ n, (3.14a) Ψ =0,whileΨ andΨ areundetermined. n n 3 1 4 | | | | n n arg(1+λ ) arg( λ )= π. (3.14b) n n 3. If n = m then = n + 1 and = m + 1, so − n+m − n+m A B N N Ψ = Ψ and Ψ = Ψ . | 1| | 4| | 2| | 3| Thesolutionλnofthisequationandhencethedecaylengthδ donotdependonthelength T oftheperiod,butonlyonthe Ineachcasetheboundaryconditionisofthecanonicalform | | ration/(n+m)=(1 √3tanϕ)/2,whichisafunctionof Ψ=(ν τ) (n σ)Ψwith − · ⊗ · theangleϕbetweenT andthearmchairorientation[seeEq. 1. ν = zˆ, n = zˆ if n > m (zigzag-type boundary (3.3)]. Inthecasen >mwhensublatticeB hasmoremodes condit−ion). thanconditions,thelargestdecaylengthδfollowsuponinter- changingnandm. 2. ν = zˆ,n = zîfn < m(zigzag-typeboundarycondi- AsseenfromFig.2, theresultingdistanceδ withinwhich tion). the zigzag-type boundary condition breaks down is zero for thezigzagorientation(ϕ = π/6)andtendsto infinityasthe 3. ν zˆ=0,n zˆ=0ifn=m(armchair-typeboundary orientationoftheboundaryapproachesthearmchairorienta- co·ndition). · tion (ϕ = 0). (For finite periodsthe divergenceis cut off at δ T a.)Theincreaseofδnearthearmchairorientation ∼| |≫ Weconcludethattheboundaryconditionisofzigzag-typefor isratherslow: Forϕ & 0.1thezigzag-typeboundarycondi- anyorientationT oftheboundary,unlessT isparalleltothe tionremainspreciseonthescaleofafewunitcellsawayfrom bonds[sothatn=mandϕ=0(moduloπ/3)]. theboundary. 5 Thedensityofedgestatesismaximalρ = 1/3aforaperfect zigzagedgeanditdecreasescontinuouslywhentheboundary orientation ϕ approaches the armchair one. Eq. (3.16) ex- plainsthenumericaldataofRef.[11],providingananalytical formulaforthedensityofedgestates. IV. STAGGEREDBOUNDARYPOTENTIAL The electron-hole symmetry (3.2), which restricts the boundary condition to being either of zigzag-type or of armchair-type,isbrokenbyanelectrostaticpotential.Herewe consider,motivatedbyRef. [3], the effectofastaggeredpo- tentialatthezigzagboundary.Weshowthattheeffectofthis potentialistochangetheboundaryconditioninacontinuous wayfromΨ= τ σ ΨtoΨ= τ (σ [zˆ n ])Ψ.The z z z B ± ⊗ ± ⊗ · × first boundary condition is of zigzag-type, while the second boundary condition is produced by an infinitely large mass termattheboundary[8]. FIG.2: Dependenceontheorientationϕofthedistanceδfromthe The staggered potential consists of a potential V = +µ, boundary within which the zigzag-type boundary condition breaks A down. Thecurveiscalculatedfromformula(3.14)validinthelimit VB = µ on the A-sites and B-sites in a total of 2N rows − |T|≫aoflargeperiods. Theboundaryconditionbecomesprecise closest to the zigzag edge parallel to the y-axis (see Fig. 3). uponapproachingthezigzagorientationϕ=π/6. Since this potential does not mix the valleys, the boundary conditionnearazigzagedgewithstaggeredpotentialhasthe form Although the presented derivation is only valid for peri- odic boundaries and low energies, such that the wavelength Ψ= τ (σ cosθ+σ sinθ)Ψ, (4.1) z z y − ⊗ ismuchlargerthanthelength T oftheboundaryperiod,we arguethattheseconditionsma|yb|erelaxed. Indeed,sincethe inaccordwiththegeneralboundarycondition(2.10). Forθ = boundaryconditionislocal,itcannotdependonthestructure 0,πwehavethezigzagboundaryconditionandforθ = π/2 ± oftheboundaryfaraway,hencetheperiodicityofthebound- wehavetheinfinite-massboundarycondition. ary cannotinfluence the boundarycondition. It can also not To calculate the angle θ we substitute Eq. (3.10) into the dependonthewavelengthoncethewavelengthislargerthan tight-binding equation (3.1) (including the staggered poten- thetypicalsizeofaboundaryfeature(ratherthanthelengthof tialattheleft-handside)andsearchforasolutioninthelimit theperiod).Sinceformostboundariesbothδandthescaleof ε = 0. Theboundaryconditionispreciseforthezigzagori- theboundaryroughnessareoftheorderofseveralunitcells, entation, so we may set αp = α′p = 0. It is sufficient to weconcludethatthezigzagboundaryconditionisingeneral consider a single valley, so we also set Ψ3 = Ψ4 = 0. The agoodapproximation. remaining nonzero components are Ψ1eiK·r ψA(i)eiKy andΨ eiK·r ψ (i)eiKy,whereiintheargu≡mentofψ 2 B A,B ≡ numbersthe unit cell away from the edge and we have used E. Densityofedgestatesnearazigzag-likeboundary that K points in the y-direction. The resulting difference equationsare A zigzag boundaryis known to support a band of disper- µψ (i)=t[ψ (i) ψ (i 1)], i=1,2,...N, (4.2a) A B B sionless states [11], which are localized within several unit − − − µψ (i)=t[ψ (i) ψ (i+1)], i=0,1,2,...N 1, cellsneartheboundary. Wecalculatethe1Ddensityofthese B A A − − (4.2b) edgestatesnearanarbitraryzigzag-likeboundary. Againas- suming that the sublattice A hasmore boundarymodesthan ψ (0)=0. (4.2c) A conditions(n < m), for each k there are (k) N lin- A A earlyindependentstates(3.8),satisfyingthNebound−arycondi- For the Ψ1,Ψ2 components of the Dirac spinor Ψ the boundarycondition(4.1)isequivalentto tion. For k = 0 the number of boundary modes is equal to 6 =n (m n)/3,sothatforeachkthereare NA − − ψA(N)/ψB(N)= tan(θ/2). (4.3) − N = (k) n=(m n)/3 (3.15) states NA − − SubstitutingthesolutionofEq.(4.2)intoEq.(4.3)gives edgestates. Thenumberoftheedgestatesforthecasewhen 1+sinh(κ)sinh(κ+2Nµ/t) n > magainfollowsuponinterchangingnandm. Theden- cosθ = , (4.4) sityρofedgestatesperunitlengthisgivenby cosh(κ)cosh(κ+2Nµ/t) Nstates m n 2 withsinhκ=µ/2t.Eq.(4.4)isexactforN 1,butitisac- ρ= = | − | = sinϕ. (3.16) ≫ T 3a√n2+nm+m2 3a| | curatewithin2%foranyN.Thedependenceoftheparameter | | 6 θoftheboundaryconditiononthestaggeredpotentialstrength Inthissectionweconsiderthemostgeneralboundarycon- µ is shown in Fig. 4 forvariousvaluesof N. The boundary dition(2.10),constrainedonlybytime-reversalsymmetry.We conditionisclosesttotheinfinitemassforµ/t 1/N,while donotrequirethattheboundaryispurelyaterminationofthe ∼ theregimesµ/t 1/N orµ/t 1correspondtoa zigzag lattice,butallowforarbitrarylocalelectricfieldsandstrained ≪ ≫ boundarycondition. bonds.TheconclusionofSec.III,thattheboundarycondition iseitherzigzag-likeorarmchair-like,doesnotapplytherefore totheanalysisgiveninthissection. The general solution of the Dirac equation (2.1) in the nanoribbonhastheformΨ(x,y)=Ψ (x)eiky. Weimpose n,k thegeneralboundarycondition(2.10), Ψ(0,y)=(ν τ) (n σ)Ψ(0,y), (5.1a) 1 1 · ⊗ · Ψ(W,y)=(ν τ) (n σ)Ψ(W,y), (5.1b) 2 2 · ⊗ · withthree-dimensionalunitvectorsν , n , restrictedbyn i i i · xˆ = 0 (i = 1,2). (Thereisnorestrictiononthe ν .) Valley i isotropyoftheDiracHamiltonian(2.2)impliesthatthespec- FIG.3: ZigzagboundarywithV = +µontheA-sites(filleddots) trumdoesnotdependonν andν separatelybutonlyonthe 1 2 and V = −µ on the B-sites (empty dots). The staggered poten- angleγ between them. Thespectrumdepends, therefore,on tialextendsover2N rowsofatomsnearesttothezigzagedge. The threeparameters: Theangleγ andtheanglesθ ,θ between 1 2 integericountsthenumberofunitcellsawayfromtheedge. thez-axisandthevectorsn ,n . 1 2 The Dirac equation HΨ = εΨ has two plane wave solu- tionsΨ exp(iky+iqx)foragivenεandk,corresponding ∝ tothetwo(realorimaginary)transversewavenumbersqthat solve (~v)2(k2 +q2) = ε2. Each of these two plane waves has a twofold valley degeneracy, so there are four indepen- dentsolutionsintotal. Sincethewavefunctioninaribbonisa linearcombinationofthesefourwaves,andsinceeachofthe Eqs.(5.1a,5.1b)hasatwo-dimensionalkernel,theseequations providefourlinearlyindependentequationstodeterminefour unknowns.TheconditionthatEq.(5.1)hasnonzerosolutions gives an implicit equation for the dispersion relation of the nanoribbon: cosθ cosθ (cosω cos2Ω)+cosωsinθ sinθ sin2Ω 1 2 1 2 − sinΩ[sinΩcosγ+sinωsin(θ θ )]=0, (5.2) 1 2 − − whereω2 =4W2[(ε/~v)2 k2]andcosΩ=~vk/ε. − For θ = θ = 0 and γ = π Eq. (5.2) reproduces the 1 2 transcendental equation of Ref. 6 for the dispersion relation ofazigzagribbon.Inthecaseθ =θ =π/2ofanarmchair- 1 2 FIG.4: Plotoftheparameterθintheboundarycondition(4.1)ata likenanoribbon,Eq.(5.2)simplifiesto zigzagedgewiththestaggeredpotentialofFig.3.Thecurvesarecal- culatedfromEq.(4.4). Thevaluesθ = 0andθ = π/2correspond, cosω =cosγ. (5.3) respectively,tothezigzagandinfinite-massboundaryconditions. This is the only case when the transverse wave function Ψ (x) is independent of the longitudinal wave number k. n,k InFig.5weplotthedispersionrelationsforseveraldifferent V. DISPERSIONRELATIONOFANANORIBBON boundaryconditions. The low energy modes of a nanoribbon with ε < ~v k | | | | Agraphenenanoribbonisacarbonmonolayerconfinedto [seepanelsa-dofFig.5]haveimaginarytransversemomen- alongandnarrowstrip. Theenergyspectrumε (k)ofthen- tum since q2 = (ε/~v)2 k2 < 0. If q becomes larger n − | | thtransversemodeisafunctionofthewavenumberk along than the ribbon width W, the corresponding wave function thestrip. Thisdispersionrelationisnonlinearbecauseofthe becomeslocalizedattheedgesofthenanoribbonanddecays confinement,which also may openup a gapin the spectrum inthebulk.Thedispersionrelation(5.2)forsuchanedgestate aroundzeroenergy. We calculatethe dependenceofthe dis- simplifiestoε=~v k sinθ forthestatelocalizednearx=0 1 persionrelationonthe boundaryconditionsatthetwo edges and ε = ~v k sin|θ| for the state localized near x = W. 2 − | | x=0andx=W ofthenanoribbon(takenalongthey-axis). These dispersive edge states with velocity vsinθ generalize 7 opposite zigzag edges have the same staggered potential, so thattheboundaryconditionis Ψ(0,y)=+τ (σ cosθ+σ sinθ)Ψ(0,y), (5.4a) z z y ⊗ Ψ(W,y)= τ (σ cosθ+σ sinθ)Ψ(W,y). (5.4b) z z y − ⊗ Thedependenceofθontheparametersµ,N ofthestaggered potentialisgivenbyEq.(4.4). Thisboundaryconditioncor- respondstoγ = π,θ = θ = θ, sothatithasagapforany 1 2 nonzeroθ. AsshowninFig.6,∆(θ)increasesmonotonically with θ from the zigzag limit ∆(0) = 0 to the infinite-mass limit∆(π/2)=π~v/W. 3 ] W2 / v h¯ [ ∆ 1 0 0 π/4 π/2 θ FIG. 6: Dependence of the band gap ∆ on the parameter θ in the staggeredpotentialboundarycondition(5.4). VI. BANDGAPOFATERMINATEDHONEYCOMB LATTICE Inthissectionwe returntothecase ofaboundaryformed purelybyterminationofthelattice.Ananoribbonwithzigzag boundaryconditionhaszerobandgapaccordingtotheDirac equation(Fig. 5a). Accordingto the tight-bindingequations FIG.5: Dispersionrelationofnanoribbonswithdifferentboundary conditions. Thelarge-wavenumberasymptotes|ε| = ~v|k|ofbulk thereis a nonzerogap∆, which howevervanishesexponen- statesareshownbydashedlines. Modesthatdonotapproachthese tially with increasing width W of the nanoribbon. We esti- asymptotesareedgestateswithdispersion|ε| = ~v|ksinθi|. The matethedecayrateof∆(W)asfollows. zigzag ribbon with γ = π and θ1 = θ2 = 0 (a) exhibits disper- Thelow energystates in a zigzag-typenanoribbonare the sionlessedgestatesatzeroenergy[11]. Ifθ1 orθ2 arenonzero(b, hybridizedzeroenergyedgestatesattheoppositeboundaries. c)theedgestatesacquirelineardispersionandifsinθ1sinθ2 > 0 Theenergyεofsuchstatesmaybeestimatedfromtheoverlap (c) aband gap opens. If γ isunequal to0or π (d) thevalleys are betweenthe edgestates localized atthe oppositeedges, ε = mixedwhichmakesallthelevelcrossingsavoidedandopensaband (~v/W)exp( W/l ). In a perfectzigzag ribbon there gap. Armchair-likeribbonswithθ1 = θ2 = π/2(e,f)aretheonly a±reedgestatesw−ithl decay=0(andε= 0),sothatthereisno ribbonshavingnoedgestates. decay bandgap. Fora ribbonwith a morecomplicatededgeshape the decay length of an edge state is limited by δ, the length within which the boundarycondition breaks down (see Sec. theknown[11]dispersionlessedgestatesata zigzagbound- III.D). This length scale provides the analytical estimate of ary(withsinθ =0). thebandgapinazigzag-likeribbon: Inspectionofthedispersionrelation(5.2)givesthefollow- ing condition for the presence of a gap in the spectrum of theDirac equationwitharbitraryboundarycondition: Either ∆ ~ve−W/δ, (6.1) the valleys should be mixed (γ = 0,π) or the edge states ∼ W 6 atoppositeboundariesshouldhaveenergiesofoppositesign withδgivenbyEqs.(3.13)and(3.14). (sinθ sinθ >0forγ =πorsinθ sinθ <0forγ =0). 1 2 1 2 Thebandgapofanarmchair-likeribbonis Asanexample,wecalculatethebandgapforthestaggered potentialboundaryconditionofSec.IV. We assumethatthe ∆=(~v/W)arccos(cosγ) (6.2) 8 2 [see Eq. (5.3) and panels e,f of Fig. 5]. Addinganotherrow ofatomsincreasesthenanoribbonwidthbyonehalfofaunit cellandincreasesγ byK R = 4π/3,sotheproduct∆W 3 · insucharibbonisanoscillatoryfunctionofW withaperiod of1.5unitcells. To test these analytical estimates, we have calculated ∆(W) numerically for various orientations and configura- ) ϕ tions of boundaries. As seen from Fig. 7, in ribbons with 1 ( f a non-armchair boundary the gap decays exponentially ∝ exp[ f(ϕ)W/a] as a functionof W. Nanoribbonswith the − sameorientationϕbutdifferentperiod T havethesamede- | | cay rate f. As seen in Fig. 8, the decay rate obtained numerically agrees well with the analytical estimate f = a/δ following from Eq. (6.1) (with δ given as a function of ϕ in Fig. 2). The numerical results of Fig. 7 are consistent with earlierstudiesofthe orientationdependenceofthebandgap 0 π/6 innanoribbons[7],buttheexponentialdecreaseofthegapfor ϕ non-armchairribbonswasnotnoticedinthosestudies. For completenesswe show in Fig. 9 ournumericalresults FIG. 8: Dependence of the gap decay rate on the orientation ϕ of forthebandgapinanarmchair-likenanoribbon(ϕ=0). We the boundary (defined in the inset of Fig. 2). The dots are the fits see that the gap oscillates with a period of 1.5 unit cells, in tonumericalresultsofthetight-bindingequations,thesolidcurveis agreementwithEq.(6.2). theanalyticalestimate(6.1). FIG.7:Dependenceofthebandgap∆ofzigzag-likenanoribbonson thewidthW.Thecurvesintheleftpanelarecalculatednumerically FIG.9: Dependenceofthebandgap∆onthewidthW foranarm- fromthetight-bindingequations.Therightpanelshowsthestructure chair ribbon (dashed line) and for a ribbon witha boundary of the oftheboundary,repeatedperiodicallyalongbothedges. sameorientationbutwithalargerperiod(solidline).Thecurvesare calculatednumericallyfromthetight-bindingequations. VII. CONCLUSION ∆ (~v/W)exp( W/δ)inananoribbonofwidthW. We ≈ − havetestedouranalyticalresultsfor∆withthenumericalso- In summary, we have demonstrated that the zigzag-type lutionofthetight-bindingequationsandfindgoodagreement. boundary condition Ψ = τ σ Ψ applies generically While the lattice termination by itself can only produce z z ± ⊗ to a terminated honeycomb lattice. The boundary condition zigzag or armchair-type boundary conditions, other types of switchesfromtheplus-signtotheminus-signatthearmchair boundaryconditionscanbereachedbybreakingtheelectron- orientationϕ=0(moduloπ/3),whentheboundaryisparal- holesymmetryofthe tight-bindingequations. We havecon- lelto1/3ofallthecarbon-carbonbonds(seeFig.10). sideredtheeffectofastaggeredpotentialatazigzagboundary The distance δ from the edge within which the boundary (producedforexamplebyedgemagnetization[3]), andhave conditionbreaksdownis minimal(= 0) atthe zigzagorien- calculatedthe correspondingboundarycondition. Itinterpo- tation ϕ = π/6 (moduloπ/3) and maximal at the armchair latessmoothlybetweenthezigzagandinfinite-massboundary orientation.Thisisthelengthscalethatgovernsthebandgap conditions,openingupagapinthespectrumthatdependson 9 thestrengthandrangeofthestaggeredpotential. Using only the Hermiticity of M, we have the 16-parameter We have calculated the dispersion relation for arbitrary representation boundaryconditionsandfoundthattheedgestateswhichare 3 dispersionlessatazigzagedgeacquireadispersionformore M = (τ σ )c , (A2) general boundary conditions. Such propagating edge states i⊗ j ij exist, for example, near a zigzagedge with staggeredpoten- iX,j=0 tial. with real coefficients c . Anticommutationwith the current ij Ourdiscoverythatthezigzagboundaryconditionisgeneric operatorbringsthisdowntothe8-parameterform explainsthefindingsofseveralcomputersimulations[11,13, 14] in which behavior characteristic of a zigzag edge was 3 observed at non-zigzag orientations. It also implies that the M = τi (ni σ), (A3) ⊗ · mechanismofgapopeningatazigzagedgeofRef.[3](pro- Xi=0 duction of a staggered potential by magnetization) applies where the n ’s are three-dimensional vectors orthogonal to i generically to any ϕ = 0. This may explain why the band n . The absence of off-diagonalterms in M2 requires that 6 B gap measurementsof Ref. [10] producedresults that did not thevectorsn , n , n aremultiplesofaunitvectorn˜ which 1 2 3 dependonthecrystallographicorientationofthenanoribbon. isorthogonalton . ThematrixM maynowberewrittenas 0 M =τ (n σ)+(ν˜ τ) (n˜ σ). (A4) 0 0 ⊗ · · ⊗ · TheequalityM2 =1furtherdemandsn2+ν˜2 =1,leading 0 tothe4-parameterrepresentation(2.8)afterredefinitionofthe vectors. APPENDIXB:DERIVATIONOFTHEBOUNDARYMODES We derive the characteristic equation (3.7) from the tight- FIG.10:Thesetwographeneflakes(orquantumdots)bothhavethe samezigzag-typeboundarycondition: Ψ = ±τz ⊗σzΨ. Thesign binding equation (3.1) and the definitions of the boundary switchesbetween+and−whenthetangenttotheboundaryhasan modes(3.4) and(3.5). In thelow energylimitε/t a/T ≪ | | anglewiththex-axiswhichisamultipleof60◦. wemaysetε 0inEq.(3.1),soitsplitsintotwodecoupled → sets of equationsfor the wave function on sublattices A and B: ψ (r)+ψ (r R )+ψ (r R )=0, (B1a) B B 1 B 2 Acknowledgments − − ψ (r)+ψ (r+R )+ψ (r+R )=0. (B1b) A A 1 A 2 ThisresearchwassupportedbytheDutchScienceFounda- SubstitutingR byR +R intheseequationsandusingthe 1 2 3 tionNWO/FOM.WeacknowledgehelpfuldiscussionswithI. definition(3.5)ofλweexpressψ(r+R )throughψ(r), 2 Adagideli,J.H.Bardarson,Ya. B.Bazaliy,andI.Snyman. ψ (r+R )= (1+λ)−1ψ (r), (B2a) B 2 B − ψ (r+R )= (1+λ)ψ (r). (B2b) A 2 A − APPENDIXA:DERIVATIONOFTHEGENERAL BOUNDARYCONDITION(2.8) Eqs.(3.5)and(B2)togetherallowtofindtheboundarymode withagivenvalueofλonthewholelattice: We first show that the anticommutation relation (2.7) fol- ψ (r+pR +qR )=λq( 1 λ)−pψ (r), (B3a) lows from the current conservation requirement (2.6). The B 2 3 − − B ψ (r+pR +qR )=λq( 1 λ)pψ (r), (B3b) current operator in the basis of eigenvectors of M has the A 2 3 A − − blockform withpandqarbitraryintegers.Substitutingψ(r+T)intoEq. (3.4) from Eq. (B3) and using T = (n+m)R +nR we X Y 1 0 2 3 nB ·J =(cid:18)Y† Z(cid:19), M =(cid:18)0 1(cid:19). (A1) arriveatthecharacteristicequation(3.7). − We now find the rootsof the Eq. (3.7) for a givenk. It is The HermitiansubblockX acts in the two-dimensionalsub- sufficientto analyze the equationfor sublattice A only since spaceofeigenvectorsofM witheigenvalue1. Toensurethat the calculationfor sublattice B is the same after interchang- Ψn J Ψ = 0foranyΨinthissubspaceitisnecessary ing n and m. The analysis of Eq. (3.7a) simplifies in polar B hand| suffi·cie|ntithat X = 0. The identity (n J)2 = 1 is coordinates, B equivalenttoYY† =1andZ =0,hence M,n· J =0. { B· } 1+λm+n = λn (B4) We now show that the most general 4 4 matrix M that | | | | × (m+n)arg( 1 λ) k narg(λ)=2πl, (B5) satisfies Eqs. (2.5) and (2.7) has the 4-parameterform (2.8). − − − − 10 with l = 0, 1, 2.... The curve defined by Eq. (B4) is tween λ and λ∗, the incrementof the left-hand side of Eq. a contour on±the±complex plane around the point λ = 1 (B5)betnweenλ∗nandλ mustbeequalto2π(n 1) 2πn − n n − ≈ which crosses points λ = 1/2 i√3/2 (see Fig. 11). (for T a),whichimmediatelyleadstoEq.(3.14)forλ . ± n − ± | |≫ Theleft-handsideofEq.(B5)isamonotonicfunctionofthe positiononthiscontour.Ifitincreasesby2π∆lontheinterval Im(λ) betweentworootsoftheequation,thenthereare∆l 1roots − inside this interval. For k = 0 both λ and λ are rootsof − + thecharacteristicequation. Sointhiscasethenumber of A N roots lying inside the unit circle can be calculated from the incrementof the left-hand side of Eq. (B5) between λ and − λ+: Re(λ) 1 0 − 1 2π 2π n m = (n+m) +n 1=n − 1. A N 2π (cid:20) 3 3 (cid:21)− − 3 − (B6) Similarly,onsublatticeB,wehave(uponinterchangingnand m), FIG.11:Plotofthesolutionsofthecharacteristicequations(B4,B5) m n =m − 1. (B7) forn = 5, m = 11, andk = 0. Thedotsaretheroots, thesolid B N − 3 − curveisthecontourdescribedbyEq.(B4),andthedashedcirclesare unitcircleswithcentersat0and−1. The same method can be applied to calculate λ . Since n therearen 1 rootsonthecontourdefinedbyEq.(B4)be- − [1] P.R.Wallace,Phys.Rev.71,622(1947). arXiv:cond-mat/0701599. [2] D.P.DiVincenzoandE.J.Mele,Phys.Rev.B29,1685(1984). [10] M.Y.Han,B.Oezyilmaz, Y.Zhang,andPh.Kim,Phys.Rev. [3] Y.-W.Son,M.L.Cohen,andS.G.Louie,Phys.Rev.Lett.97, Lett.98,206805(2007). 216803(2006). [11] K.Nakada,M.Fujita,G.Dresselhaus,andM.S.Dresselhaus, [4] E.McCannandV.I.Fal’ko,J.Phys.Condens.Matter16,2371 Phys.Rev.B54,17954(1996). (2004). [12] ThemethoddescribedinSec.IIIcanbegeneralizedtobound- [5] A.R.AkhmerovandC.W.J.Beenakker, Phys.Rev.Lett.98, arieswithN′ >n+msuchasthe“stronglydisorderedzigzag 157003(2007). boundary” of I. Martin and Ya. M. Blanter, arXiv:0705.0532. [6] L.BreyandH.A.Fertig,Phys.Rev.B73,235411(2006). Forthesenon-minimalboundariesthezigzagboundarycondi- [7] M.Ezawa,Phys.Rev.B73,045432(2006);PhysicaStatusSo- tionisstillgeneric. lidi(c)4,489(2007). [13] A.RycerzandC.W.J.Beenakker,arXiv:0709.3397. [8] M.V.BerryandR.J.Mondragon,Proc.R.Soc.LondonA412, [14] A.Rycerz,arXiv:0710.2859. 53(1987). [9] Z. Chen, Y.-M. Lin, M. J. Rooks, and Ph. Avouris,