Orbital magnetic susceptibility of finite-sized graphene

0.67 MB·English

by  Yuya Ominato

Orbital magnetic susceptibility of finite-sized graphene
Yuya Ominato and Mikito Koshino
Department of Physics, Tohoku University, Sendai 980-8578, Japan
(Dated: February 1, 2012)
Westudytheorbitalmagnetismofgrapheneribbonintheeffective-massapproximation,tofigure
out the finite-size effect on the singular susceptibility known in the bulk limit. We find that the
susceptibility at T = 0 oscillates between diamagnetism and paramagnetism as a function of εF,
in accordance with the subband structure formed by quantum confinement. In increasing T, the
2
oscillation rapidly disappears once the thermal broadening energy exceeds the subband spacing,
1
0 andthesusceptibilityχ(εF)approachesthebulklimiti.e.,athermallybroadeneddiamagneticpeak
2 centered at εF = 0. The electric current supporting the diamagnetism is found to flow near the
edge with a depth ∼ ~v/(2πkBT), with v being the band velocity, while at T = 0 the current
n distribution spreads entirely in the sample reflecting the absence of the characteristic wavelength
a in graphene. The result is applied to estimate the three-dimensional random-stacked multilayer
J
graphene,whereweshowthattheexternalmagneticfieldissignificantlyscreenedinsidethesample

1 in low temperatures, in a much stronger manner than in graphite.
3

] I. INTRODUCTION values in increasing ε , in accordance with the detailed
l F
l subband structure. In increasing temperature, the os-
a
h Graphite is known as one of the strongest diamag- cillationrapidly disappearsonce the thermal broadening
- netic materials among natural substances.1–3 This prop- energy exceeds the subband spacing, and the suscepti-
s
bility approaches bulk limit, i.e., a thermally-broadened
e erty is due to the large orbital diamagnetism related to
m the small effective mass in the band structure, i.e., nar- diamagnetic peak centered at εF = 0, independently of
the atomic configuration at the edge. We also apply a
. row energy gap between conduction and valence bands.
t similar analysis to the carbon nanotube, and find a sim-
a Thediamagneticeffectbecomesevengreateringraphene
m monolayer4–6 which is truly a zero-gap system.1,7–9 The ilar oscillation in the susceptibility.
The electric current supporting the diamagnetism
- magnetic susceptibility of graphene at zero temperature
d contains a singularity expressed as a delta function in spreads entirely in the sample at T = 0, reflecting the
n absence of the characteristic wavelength. In increasing
Fermi energyε , which divergesat Dirac point (ε =0)
o wherethetwobFandsstick,andvanishesotherwise.F1,10–17 temperature, however, the current density tends to lo-
c calize near the boundary with a depth ∼ ~v/(2πk T),
[ The orbital diamagnetism has been studied for other B
graphene-related materials, such as graphite intercala- forming an edge current circulation.
1 tion compounds,18,20–22 carbon nanotube,23–25 few-layer The analysis of the spatial distribution of the diamag-
v graphenes,19,26,27 andanorganicmaterialhavingsimilar neticcurrentisusefulinstudyingagraphenestackwhere
7
6 gapless spectrum.28 the diamagnetic current of one layer influences the elec-
tron motion in other layers. If we take a randomly-
4 In this paper, we investigate the orbitaldiamagnetism
stackedgraphenemultilayer,inwhichthe interlayercou-
6 of a graphene strip with finite width.29–41 The purpose
1. ofthisworkistwo-fold: (i)Tounderstandhowthedelta- plingis expectedto be small,45–51 the self-consistentcal-
culation shows that the diamagnetism is much stronger
0 function singularityof the bulk limit is relaxedin a real-
thaningraphite,andtheexternalmagneticfieldissignif-
2 istic finite-sized graphene system. In the literatures, the
1 orbital susceptibility of graphene nanoribbons was cal- icantly screened inside the sample in low temperatures.
v: culated for the Fermi energies near the Dirac point,31,44 Thepaperisorganizedasfollows. InSec.II,webriefly
revieweffectivemassdescriptionofelectronsingraphene
i whilethebehavioroffDiracpointandtherelationtothe
X ribbon,andformulatetheorbitalmagneticsusceptibility.
bulk susceptibility is not well understood. (ii) To study
r the diamagnetic current flow on graphene. In the con- We presentthe numericalresults anddetailed discussion
a
in Sec. III, as well as a similar analysis for the carbon
ventionaldiamagnetismofmetal, we usually expectthat
nanotube in Sec. IV. We argue the diamagnetism of
the current circulates near the surface with a depth of
randomly stacked graphene multilayers in Sec. V. The
the order of the Fermi wave length λ . In graphene, the
F
conclusion is given in Sec. VI.
only characteristic length scale λ intrinsically diverges,
F
and we expect a peculiar current distribution different
from the conventional system.
II. FORMULATIONS
To address above problems, here we calculate the
orbital susceptibility and the current distribution of
graphene ribbon with an arbitrary width, in various A. Effective mass approximation
Fermienergiesε andtemperaturesT,usingtheeffective
F
mass approximation. We find that the susceptibility at Grapheneiscomposedofahoneycombnetworkofcar-
T = 0 oscillates between diamagnetic and paramagnetic bonatoms,whereaunitcellcontainsapairofsublattices,
2
(a) Zigzag ribbon wave amplitudes are written as53
y
ψ (R )=eiK·RAFK(R )+eiηeiK′·RAFK′(R ),
A B A A A A A A
ψ (R )=−ωeiηeiK·RBFK(R )+eiK′·RBFK′(R ).
B B B B B B
(2)
in terms of the slowly-varying envelope functions
A B FK,FK,FK′,and FK′. The envelope functions satisfy
b A B A B
x the Schr¨odinger equation,1,7–9,52,53
η
a H F(r)=εF(r), (3)
0
x=0 x=Lx with
0 kˆ −ikˆ 0 0
(b) Armchair ribbon x y
kˆ +ikˆ 0 0 0
y H0 =~v x 0 y 0 0 kˆ +ikˆ , (4)
x y
 0 0 kˆ −ikˆ 0 
 x y 
A A  
FK(r)
A
B B FK(r)
b a x F(r)=FBK′(r), (5)
A A
FK′(r)
 B 
 
where kˆ =−i∇ and v is the band velocity.
B
The electronic states of the graphene ribbon can be
x=0 x=Lx correctly described by setting the appropriate boundary
condition to the effective mass Hamiltonian.34 Now the
eigenstates are labeled by k since the system is transla-
y
FIG.1: Atomicstructuresofgrapheneribbonswith(a)zigzag tionally symmetric along y. A wave function of ky and
boundary and (b) armchair boundary, respectively. Dashed the energy ε is generally written as
circles indicate missing sites beyondthe boundary.
Aeikxx+Be−ikxx
F(r)=eikyys Aei(kxx+θ)−Be−i(kxx+θ) , (6)
Ceikxx+De−ikxx
denoted by A and B. Fig.1 (a) and (b) show the atomic (cid:0) (cid:1)
s Cei(kxx−θ)−De−i(kxx−θ) 
structure of zigzag and armchair graphene ribbons, re-  
spectively, where a and b are primitive translation vec-  
(cid:0) (cid:1)
torsofinfinitegraphene. Thelatticeconstantisgivenby wherek2 =ε2/~2v2−k2,eiθ =(k +ik )/ k2+k2,s=
x y x y x y
a = |a| ≈ 0.246 nm. For both cases we set y-axis along
ε/|ε|, and A,B,C and D are numbers toqbe determined
the ribbon, and set x = 0 and L to the line of missing
x bysatisfying the boundarycondition,as wewillarguein
sites nearest from the edge. We define η as the angle
the following.
between x axis and a, which is π/6 for zigzag, and 0 for
armchair boundary.
In a tight-binding model, the wave function of B. Zigzag boundary
graphene electron is written as
In the zigzag ribbon, the boundary condition is given
ψ(r)= ψA(RA)φ(r−RA)+ ψB(RB)φ(r−RB), by ψA(RA) = 0 at x = 0, and ψB(RB) = 0 at x = Lx.
By using Eq. (2), this is translated to the condition for
XRA XRB (1) the envelope function as
FK(0,y)=0,
where R and R are the positions of A-sites and B- A
A B FK(L ,y)=0,
sites,respectively,andφ(r)denotesthe wavefunctionof B x
the p orbital of a carbon atom. FK′(0,y)=0,
z A
For states in the vicinity of the Fermi level ε =0, the FK′(L ,y)=0, (7)
B x
3
15 The wavefunction and the energy then becomes
k2
10 ttaannh k xkLxxLx y (=k ykLxx/ k< y 1) eikyy −ssiinshinκh(κx0x−L )
k1 Fs0ky(r)=A1 L L  00 x ,
x y
y 5 k0 κ0 k1 y( k=y Lkxx />k y1) k2 A0 = −1+ sipnh2κ0Lx −1/2, 0 
2κ L
0 (cid:18) 0 x (cid:19)
ε =s~v −κ2+k2. (12)
s0ky 0 y
-5 q
0 π/Lx 2π/Lx This actually describes the edge state localized at the
kx boundary x = 0 and L giving a nearly flat energy
x
band.29–31
The eigenenergy ε represents the n-th branch of
FIG. 2: Plots to findsolutions of Eq.(9) and Eq.(11). conduction (s = +)snaknyd valence (s = −) bands re-
spectively. Energy band structure of K as a function
ofk isshownassolidcurvesinFig.3(a). Eigenstatesfor
y
whichkeepsthestatesatKandthoseatK′ independent. K′ point are obtained similarly, where the energy band
For an eigenstate for the K point, we apply the first two structure is equivalent to Fig. 3(a) with ky inverted to
lines of Eq.(7) to Eq.(6), to obtain −ky. The flat band of edge states of K and K′ are con-
nected in a wave number away from K or K′.29–31
A+B =0,
s Aei(kxLx+θ)−Be−i(kxLx+θ) =0. (8) C. Armchair boundary
(cid:16) (cid:17)
In the armchair ribbon, the boundary condition im-
To have a solution other than A=B =0, we require34 poses both of ψ (R ) = 0 and ψ (R ) = 0 at each of
A A B B
x=0 and x=L . The corresponding conditions for the
x
k envelope functions are written as
x
k = . (9)
y
tank L
x x FK(0,y)+FK′(0,y)=0,
A A
FK(0,y)−FK′(0,y)=0,
B B
For given k , we define k (n = 0,1,2,···) as solution
of Eq.(9) in kxy satisfying nπn<knLx <(n+1)π. Corre- FAK(Lx,y)+ω−2NFAK′(Lx,y)=0,
sponding eigenstates and the energy are obtained as FK(L ,y)−ω−2NFK′(L ,y)=0, (13)
B x B x
isinknx where N = Lx/a is the number of honeycomb lattices
F (r)=A eikyy s(−1)n+1sinkn(x−Lx) , between x = 0 and Lx, which can be integer or half-
snky n L L  0  integer depending on the position of the edge.
x y
0 Applying above conditions to Eq.(6), we obtain
 
p  
sin2knLx −1/2 1 1 1 1 A 0
A = 1− ,
n (cid:18) 2knLx (cid:19)  eeiiλθ −ee−−iλiθ −αee−iλiθ αee−iθiλ BC=00,
ε =s~v k2 +k2, (10)
snky n y ei(λ+θ) −e−i(λ+θ) −αei(λ+θ) αe−i(λ+θ)D 0
q    (14)
for n = 0,1,2,··· . For the normalization of the wave-
function,weassumedtheperiodicboundaryconditionin whereα=ω−2N andλ=k L . The determinantofma-
x x
y-direction with a large enough period L . trixinEq.(14)shouldvanishtohaveanon-zerosolution.
y
Solutionk isobtainedbysearchingforcrossingpoints This condition is reduced to
n
of tankxLx and kx/ky as illustrated in Fig. 2. When π ν
kyLx >1,the firstsolutionk0 becomesapureimaginary kx =kn ≡ L n− 3 , n=0,±1,±2,··· , (15)
x
number iκ0 which satisfies (cid:16) (cid:17)
ν is an integer (0,±1) defined by
κ
0
k = . (11)
y tanhκ L 2N =3m+ν, (16)
0 x
4
with integer m. The eigenstate and energy are obtained case in armchair ribbons. This is obvious from at the
as matrixelementofˆj betweentwoeigenstatesF andF′,
y
eiknx hF′|ˆj (r)|Fi=iev F′K(r)∗FK(r)−F′K(r)∗FK(r) .
Fsnky(r)= 2 eiLkyyL  s−eie(k−nixk+nxθ) , y (cid:2) A B B A (2(cid:3)3)
x y
se−i(knx−θ) Inthe wavefunctionofzigzagribbon,Eqs.(10)and(12),
ε =s~v pk2 +k2.  (17) the component FAK is zero atx=0, andFBK is at Lx, so
snky n y that Eq. (23) vanishes at the both edges.
q The currentdensity onxy-plane is related to the local
When ν = 0, the energy bands of n = 0 and s = ± magnetic moment m(r) in z-direction by
sticktogetherandthus the systemis metallic,while oth-
erwise a gap opens at zero energy and the system be- ∂m ∂m
j =c , j =−c . (24)
comes a semiconductor. Energy bands for metallic arm- x ∂y y ∂x
chairribbon(ν =0)andsemiconductingarmchairribbon
(ν =±1)are shownin Fig.3(b) and(c), respectively. In In the present case, m(r) depends only on x so that the
(c), the labeling n is for the case of ν = +1, while n total magnetization per area is
becomes −n in ν =−1.
1 Lx Ly
M = m(x)dxdy
L L
x y Z0 Z0
D. Orbital susceptibility 1 Lx L
x
= x− j (x)dx. (25)
y
cL 2
x Z0 (cid:18) (cid:19)
To calculate the orbital diamagnetism, we consider a
graphene ribbon under a uniform magnetic field B per- The magnetic susceptibility is written as
pendiculartographeneplane. WetaketheLandaugauge
2
and set the vector potential as M 2 (δH/B)αβ
χ= lim = f(ε ) .
α
A(r)= 0, B x− Lx . (18) B→0 B cLxLy Xα βX(6=α)(cid:12)(cid:12) εβ −εα (cid:12)(cid:12)(26)
2
(cid:20) (cid:18) (cid:19)(cid:21)
We calculate Eqs. (22) and (26) numerically. As we
The Hamiltonian in presence of the magnetic field is ob-
have infinite energy bands below zero, we introduce a
tained by replacing kˆ by kˆ+eA/(~c), as
cut-offfunctiong(ε )whichsmoothlyvanishes|ε |>ε .
α α c
e In the following, we take εc =50ε0 where
H=H +δH, δH= vˆA , (19)
0 y y
c 2π~v
ε = (27)
0
where c is the light velocity, and L
x
isthetypicalenergyscaleforthesubbandstructure. The
0 −i 0 0
1 ∂H i 0 0 0 result is actually converging in the limit of large εc.
vˆy = ~∂kˆ =v0 0 0 i. (20) Thesusceptibility ofthe infinite bulk grapheneatzero
y temperature is given by1,13,20
0 0 −i 0
 
 
e2v2
Theoperatoroftheelectriccurrentdensityisgivenby χ (ε )=−g g δ(ε ), (28)
gr F v s6πc2 F
e
ˆjy(r)=−2{vˆyδ(r−r′)+δ(r−r′)vˆy}. (21) where gv = gs = 2 are the valley and spin degeneracies,
respectively. At finite temperature, it becomes
In the first order perturbation in δH, the expectation
∞ ∂f(ε)
value of the current density is written as
χ (µ;T)= dε − χ (ε)
gr gr
∂ε
Z−∞ (cid:18) (cid:19)
2Re (δH) (ˆj (r)) e2v2 1
αβ y βα
jy(r)= f(εα) h ε −ε i, (22) =−gvgs24πc2kBT cosh[µ/(2kBT)]. (29)
α β
Xα βX(6=α)
In the graphene ribbon, the characteristic unit of the
where α and β represent the unperturbed eigenstates of susceptibility can be chosen as
grapheneribbon, andf(ε)=1/[1+eβ(ε−µ)]is the Fermi
distribution function with the chemical potential µ. e2v2 1
χ =g g . (30)
Inazigzagribbon,thecurrentdensityalwaysvanishes 0 v s6πc2ε
0
at the edges x = 0 and L , while it is not generally the
x
5
(a) Zigzag ribbon (b) Armchair ribbon (ν=0) (c) Armchair ribbon (ν= +− 1)
2 s= − s=+ s= − s=+ s= − s=+
)x
L
π/
2
of 0
s
nit
u
(y 2 1 1 2
k-2 2 1 0 n=0 1 2 +−2 +−1 0 n=0 +−1 +−2 -2 -1 0 n=0 -1 -2
2
)0 1
χ
of
s
nit 0
u
χ(
-
-1
-1 0 1 -1 0 1 -1 0 1
εF (units of ε0) εF (units of ε0) εF (units of ε0)
FIG.3: Bandstructure(upperpanel)andmagneticsusceptibilityasafunctionofεF (lower),of(a)zigzag,(b)metallicarmchair
(ν = 0) and (c) semiconducting armchair (ν =±1) graphene ribbons. In upper panels, solid (black) and dashed (red) curves
indicatethebandstructuresat zeroandafinitemagnetic field,respectively. Forthelatter,theenergy bandiscalculated with
the perturbation theory in a magnetic field B. where we take B = B0 for (a), and B = 0.5B0 in (b) and (c) for illustrative
purpose, and B0 =(h/e)/L2x.
Armchair ribbon (ν=0) III. NUMERICAL RESULTS
2
T / T0 = 0 A. Magnetic susceptibility
0.16
χ)01 0.80 Lower panels of Fig.3 show the orbital susceptibility
s of χ(εF) of (a) zigzag, (b) metallic armchair (ν = 0) and
nit (c) semiconducting armchair (ν = ±1) ribbons at zero
u 0
χ( temperature, where the upward direction represents the
- negative (i.e., diamagnetic) susceptibility. The figures
are to be compared with the band structures in upper
-1 panels. In every case, the magnitude of χ becomes the
maximum at ε = 0, and oscillates as a function of ε
F F
-4 -2 0 2 4 in accordance with the subband structure. In the posi-
µ (units of ε0) tive energy region, for example, the curve sharply rises
whenasubbandstartstobeoccupiedbyelectrons,while
it tends to decrease otherwise. In large |ε |, the ampli-
F
FIG. 4: Magnetic susceptibility against chemical potential tude of the oscillation slowly attenuates approximately
µ for metallic armchair ribbon (ν = 0) at several different in proportionalto 1/ |ε |.
F
temperatures. The oscillating feature can be understood in terms of
p
the band energy shift in an infinitesimal magnetic field.
In the upper figures of Fig. 3, we plot as broken curves
theenergybandinsomesmallBcalculatedbythesecond
orderperturbation. Here the amplitude B is set to some
finitevalueforillustrativepurpose. Generallythesystem
isdiamagneticwhenthetotalenergyshiftcausedbyB is
positive,andparamagneticwhennegative. Inthe metal-
6
)x 20 Armchair ribbon (ν=0)
χcB / L0 10 ZAirgmzcahga riirb rbibobnon (ν=0) 2 20 χB/ L) x0
nitis of 0 Armchair ribbon (ν= +− 1) 1 0 nitsof c
u u
y ( )0 -20 j(y
ensit -10 εs of 0
urrent d -20 ε(unit F
C 0 0.2 0.4 0.6 0.8 1.0
-1
x coordinate (units of L )
x
-2
FIG. 5: Diamagnetic current density jy(x) of different types
0 0.2 0.4 0.6 0.8 1.0
of graphene ribbons with εF =0 at T =0. x coordinate (units of Lx)
lic armchair ribbon (b), for example, we see that a pair FIG. 6: Two-dimensional density plot jy(x;εF) of the dia-
magneticcurrentdensityofmetallicarmchairribbon(ν =0),
of the first subbands (n = 0,s = ±) shift towards zero
asafunctionofpositionx(horizontalaxis)andFermienergy
energy,duetothelevelrepulsionsfromexcitedsubbands
(vertical).
nearby. Allotherbands(n>0)moveinthe oppositedi-
rection away from zero energy, while the absolute shifts
Armchair ribbon (ν=0)
are much smaller than that of n=0.
lenWcehebnanεdF(sis=ze−ro,,nth=e 0e)neerxgcyeedgsaintheofentheregyfirlsotssvao-f B / L)x 20 T / T00 .=1 06
0
all other valence bands (s = −,|n| ≥ 1), resulting in χc 10 0.80
the total diamagnetism. When ε is shifted to positive of
F s
side, the diamagnetism decreases because the first con- niti 0
gdiuvcetsiopnarbaamnadgn(seti=sm+. ,Wnh=en0t)hehsaescoanndecgoantdivuectsihoinftbaanndd sity (u hv/(2πkBT)
n -10
(s=+,|n|=1)startsto be filled,the susceptibility sud- e
d
denly jumps to diamagnetic direction, because the shift nt
e -20
ispositivethereandalsothedensityofstatesdivergesat urr
C
the bandbottom. The oscillationofother types, (a) and
0 0.2 0.4 0.6 0.8 1.0
(c), can be explained in a similar manner. x coordinate (units of L)
x
InFig.4,thesusceptibilityatseveraldifferenttemper-
atures is plotted as a function of the chemical potential
µ. We here choose the metallic armchair ribbon (ν = 0) FIG. 7: Diamagnetic current density jy(x) of metallic arm-
while the qualitativepropertyis the samein othercases. chairribbon(ν =0)withεF =0,atseveraldifferenttemper-
atures. The vertical arrows indicate the characteristic length
We define the characteristic temperature scale T as
0 scale ~v/(2πkBT) measured from x=0.
k T =ε , (31)
B 0 0
B. Diamagnetic current density
atwhichthethermalbroadeningenergyisoftheorderof
the subband interval energy. We see that the oscillation Fig. 5 shows the diamagnetic current density j (x) in
y
rapidly disappears once T becomes of the order of T , differenttypes ofgrapheneribbons atε =0 andT =0,
0 F
leaving only single diamagnetic peak at ε = 0. When calculated in the first order perturbation of B. The unit
F
T >∼ T , the curve becomes almost identical with the ofcurrentdensityistakenascχ B/L . Thecurrentflows
0 0 x
bulk susceptibility, Eq.(29). We also confirmed that the inoppositedirectionsintheleft-handsideandright-hand
integration of χ over ε is identical with the bulk value side of the ribbon, to make a magnetization perpendicu-
F
−g g e2v2/(6πc2) within the numerical accuracy, for all lar to the layer. Reflecting the absence of the character-
v s
the types of ribbons considered here. This fact suggests isticlengthscale,thecurrentdistributionisnotlocalized
that the finite size effect disappears T >∼T , regardless to the edge but spread in the entire width in a form of
0
of the edge configuration. slowly-varying monotonic function.
7
Inzigzagribbons,thecurrentdensityactuallybecomes agnetism is given by
absolute zero at x = 0 and Lx, in accordance with the g
constraint argued in the previous section. The current χpara = 2 µ2BD(ε), (33)
sharply drops to zero at the edges, and some oscillatory (cid:16) (cid:17)
where g ∼ 2 is the g-factor for graphene electron, µ =
featureremainsaroundtheedgeduetoafinitecut-offen- B
e~/(2mc) is the Bohr-magneton with m being the free-
ergy. When we increase the energy cut-off (not shown),
electron mass, and D(ε) is the density of states per area
the curve appears to slowly approach a fixed curve hav-
given by the zero-energy flat band. Since the number of
ing a discontinuous jump at the edges. In the armchair
edge states accommodatedin a ribbon of the length L is
ribbons, j is notnecessarilyzerobutlogarithmicallydi-
y ∼L/a per spin and per valley,31 we have
vergesatthe bothedges. The numericalcalculationcon-
vergesmuchmorerapidlythere,sincethere is nodiscon- g g L
v s
D(ε)∼ δ(ε), (34)
tinuity as in the zigzag case.
L2 a
Fig. 6 is the two-dimensional density plot j (x;ε ) of
y F which gives a delta-function singularity in χ .
thediamagneticcurrentdensityofmetallicarmchairrib- para
By comparing χ with the bulk orbital diamag-
bon, as a function of position x (horizontal axis) and para
netism χ , Eq. (28), we obtain
Fermi energy (vertical). In increasing ε , the current dia
F
distribution begins to oscillate as a function of x, with a χ 3π ~2 a
para
characteristic wave length of the order of kF =εF/(~v). χ = 2 m2v2aL ∼1.0× L, (35)
dia
The temperature dependence ofthe currentdensity at (cid:12) (cid:12)
which is neg(cid:12)ligible(cid:12)in a wide strip with L ≫ a. In a low
ε = 0 is shown in Fig.7 for the same metallic arm- (cid:12) (cid:12)
F
temperaturesuchthatk T <∼ε ,however,χ cannot
chair ribbon. When T becomes as large as T , the cur- B 0 dia
0
be regarded as thermally broadened delta-function due
rent distribution is localized at the boundary forming
to the effect of the subband formation, and then χ
the counter edge currents. This is the same temperature para
overcomes χ only at ε =0.
range where the oscillation of χ disappears and the bulk dia F
limit is achieved. The depth of the current distribution
is characterized by
IV. CARBON NANOTUBES
~v
λ (T)= , (32)
edge The carbon nanotube is a quasi-one-dimensional
2πk T
B
system similar to graphene ribbon, but different in
which shrinks in increasing temperature. With the band that there are no edges.42,43 Experimentally, graphene
velocity of graphene, v ≈ 106 m/s, it is estimated as nanoribbons with smooth edges can be obtained by un-
λ ≈[1/T(K)]µm. zipping the carbon nanotubes, i.e., lengthwise cutting of
edge
This behavioris intuitively explainedusingthe plotof carbonnanotube side walls.40,41 Thenwe mayaskwhich
j (x;ε ) in Fig.6. The current density at a finite T is of the ribbon and the original nanotube has greater dia-
y F
obtained by integrating jy(x;εF) in εF with the thermal magnetism, and how the susceptibility oscillation in εF
averagingfactor−∂f/∂ε. Thecurrentofthemiddlepart changes in unzipping. The orbital susceptibility of car-
of the ribbon vanishes in averaging the oscillating func- bon nanotube was theoretically studied for small Fermi
tion in ε , while the cancellation is not complete only energies in the effective mass approximation23,24. Here
F
near the edges, since j is always positive and negative we compute full Fermi energy dependence in parallel
y
in left and right ends, respectively. The similar temper- fashion to the analysis for ribbons.
ature dependence of the current distribution is found in Acarbonnanotube is characterizedby a chiralvector,
other types of ribbons considered here. This suggests L=n a+n b, (36)
a b
that,inanyfinitepiecesofgraphenewithlengthscaleL,
the finite-size effect disappears when T > T , and then where the atom at L on a graphene sheet is rolled up
0
the diamagnetic currentcirculates only near edge with a onto the origin in constructing a tube. The bound-
depth λ . ary condition is given by ψ (R ) = ψ (R +L) and
edge A A A A
ψ (R ) = ψ (R +L). For the effective mass wave-
B B B B
function, it is written as53
C. Relation to spin paramagnetism 2πi
FK(r+L)=exp − ν FK(r),
3
(cid:18) (cid:19)
We neglect the effect of the electron spin through out
thepresentanalysis. Inazigzagribbon,particularly,the FK′(r+L)=exp +2πiν FK′(r). (37)
3
largedensityofstatescontributedbythezero-energyflat (cid:18) (cid:19)
bandis expectedto giveasignificantmagnitudeofPauli Here FK = (FK,FK) etc., and ν is an integer (0,±1)
paramagnetism and reduce the orbital diamagnetism.31 defined by A B
The ratio between two effects can be quantitatively
n +n =3m+ν, (38)
estimated as follows. The susceptibility of Pauli param- a b
8
with integer m. wherek =|ε |/(~v)istheFermiwavenumberandθ(x)
F F
ForK point,the eigenstatesareimmediately obtained is defined by θ(x) =1(x>0) and 0(x< 0). Using Eqs.
as (42) and (43), the flux-averaged susceptibility integrate
is shown to be
eiknx
FsKnky(r)= e2iLkyxyLy sei(k00nx+θ), DZ−∞∞χtube(εF)dεFEϕ = 12(cid:18)−gvgs6eπ2vc22(cid:19), (44)
p  
  which is exactly half of graphene’s, suggesting that the
ε =s~v k2+k2, (39)
snky y n susceptibility is effectively smaller in nanotube than in
q ribbon. This is simply because the B-field component
where
penetrating the lattice plane is smaller in the nanotube
2π ν
k ≡ n− , n=0,±1,±2,···, (40) due to its cylindrical shape.
n
Lx 3 The susceptibility before taking flux average can be
(cid:16) (cid:17)
and L = |L| and L is length of the carbon nanotube. calculated in numerics. Fig. 8 (a) shows χtube(εF) for
x y
the metallic (ν = 0) and the semiconducting (ν = ±1)
The system is metallic when ν =0, and semiconducting
nanotubes,togetherwiththefluxaverage. Ithasanoscil-
when ν =±1. The band structure looks similar to arm-
latorybehaviorsimilartothe grapheneribbon’s,whileχ
chairgrapheneribbon’s,buttheunitofmomentumquan-
in |ε |>~vq/2 completely vanishes after flux average.16
tization doubled compared to Eq. (15), leading to wider F
In increasing temperature (not shown), the oscillation
energy spacing between subbands. The energy band for
K′ is obtained by replacingk by −k and also ν by −ν. immediately disappears, leaving a single peak regardless
y y
of ν, similar to the graphene ribbon.
When a uniform magnetic field B is applied perpen-
Fig.8(a)comparesthe susceptibility ofa carbonnan-
dicularly to the nanotube axis, the vector potential can
otube and that of corresponding graphene ribbon un-
be taken as
zipped from the same nanotube. Here we chose a zigzag
BL 2πx
A(r)= 0, x sin . (41) ribbonasanexample,whenthe correspondingnanotube
(cid:18) 2π Lx (cid:19) becomes an armchair nanotube which is always metal-
We should note that the expression differs from that for
ribbon, Eq. (18), because the magnetic field perpendicu-
lartothe tube surfaceisnotaconstant,buta sinusoidal (a) CNT (ν=0)
function in x. Except for that, the magnetic suscepti- CNT (ν= +− 1)
)00.4
bility χ (ε ) is calculated in the same formula, Eq. χ CNT (av.)
(26). tube F of
s
The susceptibility of the carbonnanotube is naturally nit
relatedtothatofgrapheneagainstaspatialvaryingmag- u 0
(
netic field B(q)sinqx with q = 2π/L ≡ q . When χ
x 0 -
we define the q-dependent susceptibility of graphene as
χgr(q)≡m(q)/B(q),16 we obtain a relation -0.4
1
hχtube(εF)iϕ = 2χgr(q0;εF). (42) (b) CNT(ν=0)
1
) Zigzag ribbon
Here hiϕ represents an average over a phase factor ϕ χ0
whichtwiststhe boundaryconditionofcarbonnanotube of
as ψ(r+L)=exp(2πiϕ)ψ(r). Physically,the phase fac- nits
tor corresponds to threading a magnetic flux of (h/e)ϕ u 0
(
into the nanotube cross section.23,53 It changes momen- χ
-
tum quantization of Eq. (40) to k = (2π/L )(n+ϕ−
n x
ν/3),andtheflux averagingoverϕsmearsthedifference
in ν. The factor 1/2in Eq.(42) enters because the aver- -1
-2 -1 0 1 2
ageofthesquaredmagneticfieldonthenanotubesurface εF (units of ε0)
is B(q )2/2.
0
At the zerotemperature, χ (q;ε ) is explicitly evalu-
gr F
ated as16
FIG. 8: (a) Magnetic susceptibility of carbon nanotubes.
gvgse2v 1 Solid (black), dashed (red) and dotted (green) curves are for
χgr(q;εF) = − 16~c2 q θ(q−2kF) metallic (ν = 0), semiconducting (ν = ±1), and the flux av-
erage, respectively. (b) Magnetic susceptibility of the carbon
× 1+ 22kF 1− 2kF 2− 2 sin−1 2kF ,(43) nanotube of ν = 0 (solid black) and the zigzag ribbon un-
" π q s q π q # zipped from thesame nanotube (dashed red).
(cid:16) (cid:17)
9
lic (ν = 0).53 The oscillation period of the nanotube is Since the ring current I exists every interlayer distance
approximately twice as large as that of the ribbon, re- d, it induces a counter-magnetic field inside the system
flecting the wider subband spacing. Overall magnitude as
ofχissmallerinnanotuberoughlybyfactor2. Theinte-
4πI
grateofsusceptibilityinεF differsinfactor2innumerical ∆B = . (47)
accuracy, in accordance with the above arguments. c d
Solving the set of equations, we find that the dimen-
sionless volume susceptibility becomes
V. RANDOMLY-STACKED MULTILAYER
GRAPHENE ∆B −1
χ ≡ = . (48)
3D
B 1−d/(4πχ )
ext gr
Thediamagnetismcanbe madeevengreaterbystack-
Atthe chargeneutralpointµ=0,inparticular,wehave
ing graphenes in three dimensions. The recent experi-
mentaltechniquerealizesanovelkindofgraphenemulti-
−1
layer in which successive layers are stackedwith random χ (µ=0)= , (49)
3D
1+k T/∆
rotating angles.45–47 There it is known that the inter- B
layer coupling is significantly weakened and the Dirac
where ∆ is a characteristic energy scale defined by
coneiskeptalmostintactnearzeroenergyaslongasthe
rotating angle is not too small.48–51 The orbital suscep- g g v 2 e2
v s
tibility of such a system is expected to be muchstronger ∆= ≈0.03meV, (50)
6 c d
thangraphiteinwhichthedelta-functionpeakofχ(εF)is (cid:16) (cid:17)
much broadened and shortened by the regular interlayer anddisassumedtobetheinterlayerspacingofgraphite,
coupling.18,19 0.334 nm.
Here we consider the orbital diamagnetism of a finite- In decreasing the temperature, χ3D monotonically in-
sized piece of random-stacked graphene multilayers. In creases in the negative direction, and approaches −1,
calculations,we self-consistently include the effect of the where the prefect magnetic field screening is achieved.
counter magnetic field induced by the diamagnetic cur- This reflects the property of the single-layer susceptibil-
rent itself. This is essential because, as we will show in ity,Eq.(29),ofwhichpeakvalue atDiracpointdiverges
the following, the counter magnetic field of this system inT →0. Incontrast,χ3D ofthegraphiteisofthe order
can be of the same order of the external magnetic field, of10−4andisnotmuchenhancedinlowtemperatures,57
and even nearly perfect screening is possible in low tem- becauseχ(εF)isalreadybroadenedbytheinterlayercou-
peratures. plingenergyabout4000K.18,19Athree-dimensionalbulk
For simplicity, we completely neglect the interlayer materialcomposed of random-stackedgraphenes,if real-
coupling and regard the system as a set of independent ized, would be the strongest diamagnetic material than
single layer graphenes. We also assume that each layer anyother knownsubstances exceptfor the superconduc-
has the identical shape with a characteristic length scale tors. Including the effect of the residual interlayer cou-
L, and that the system is large enough that the thermal plingbetweenmisorientedlayersmaysetthe upperlimit
broadeningenergykBT ismuchlargerthan2π~v/L. Ac- to χ3D, while we leave the detailed analysis for a future
cording to the previous discussions, we then expect that problem.
the susceptibility of each layer is given by the bulk limit
χ in Eq. (29), and also the depth of the edge current
gr
λ of Eq. (32) can be neglected with respect to the VI. CONCLUSION
edge
system size L.
LetusconsiderasituationwhereaexternalfieldB is Wehavestudiedtheorbitaldiamagnetismofgraphene
ext
appliedperpendicularlytographeneplaneoftherandom ribbonsusingtheeffective-massapproximation,tofigure
stackedmultilayer. The total magnetic field B penetrat- out its dependence on temperature and Fermi energy,
ing the system is and also the finite-size effect on the delta-function sin-
gularity in the bulk limit. In increasing temperature, an
B =B +∆B, (45) oscillatory behavior in the orbital susceptibility χ(ε ) is
ext F
eventuallysmearedout,approachingthebulklimiti.e.,a
where ∆B is the counter field caused by graphene elec- thermally broadened delta-function centered at ε = 0.
F
trons. The total field B induces the magnetism in each The electric current responsible for the diamagnetism
layer, M = χ BS, with S being the area of the layer. spreads entirely in the sample at T = 0 reflecting the
gr
This is relatedto the diamagnetic edge currentI ofeach absence of the characteristic wavelength, while as T is
single layer by increased, the current density tends to localize near the
edge with a depth ∼~v/(2πk T).
B
cM We alsoseeacarbonnanotube,anotherformofquasi-
I = =cχ B. (46)
gr
S one-dimensional carbon, exhibits a similar oscillation in
10
χ(ε ), but the overall magnitude is reduced by a fac- ACKNOWLEDGMENTS
F
tor 2 compared to the corresponding ribbon having the
same width. The result is applied to estimate the three-
dimensionalbulksusceptibilityofrandom-stackedmulti-
layergraphene. There weshowedthat the externalmag- This project has been funded by JST-EPSRC Japan-
netic field is significantly screened inside the sample. UK Cooperative ProgrammeGrant No. EP/H025804/1.
1 J. W.McClure, Phys.Rev.104, 666 (1956). haus, Phys. Rev.B 54, 17954 (1996).
2 J. W.McClure, Phys.Rev.119, 606 (1960). 31 K.Wakabayashi,M.Fujita,H.Ajiki,andM.Sigrist,Phys.
3 M. P. Sharma, L. G. Johnson, and J. W. McClure, Phys. Rev.B59, 8271 (1999).
Rev.B 9, 2467 (1974). 32 M. Ezawa, Phys.Rev.B 73, 045432 (2006).
4 K.S.Novoselov, A.K.Geim,S.V.Morozov, D.Jiang, Y. 33 L. Brey, and H.Fertig, Phys.Rev. B 73, 195408 (2006).
Zhang,S.V.Dubonos,I.V.Grigorieva, andA.A.Firsov, 34 L. Brey, and H.Fertig, Phys.Rev. B 73, 235411 (2006).
Science 306, 666 (2004). 35 Y-W.Son,M.L.Cohen,S.G.Louie,Phys.Rev.Lett.97,
5 K.S.Novoselov,A.K.Geim,S.V.Morozov,D.Jiang, M. 216803 (2006).
I. Katsnelson, I. V. Grigorieva, S. V. Dubonos, and A. A. 36 Y-W. Son, M. L. Cohen, S. G. Louie, Nature 444, 347
Firsov, Nature438, 197 (2005). (2006).
6 Y.Zhang, Y.-W.Tan, H.L.Stormer,and P.Kim, Nature 37 M. Y. Han, B. O¨zyilmaz, Y. Zhang, and P. Kim, Phys.
438, 201 (2005). Rev.Lett. 98, 206805 (2007).
7 J. C. Slonczewski and P. R. Weiss, Phys. Rev. 109, 272 38 Z.Chen,Y.-M.Lin,M.J.Rooks,andP.Avouris,Physica
(1958). E 40, 228 (2007).
8 D. P. DiVincenzo and E. J. Mele, Phys. Rev. B 29, 1685 39 X.Li,X.Wang,L.Zhang,S.Lee,andH.Dai,Science319,
(1984). 1229 (2008).
9 G. W. Semenoff, Phys. Rev.Lett. 53, 2449 (1984). 40 D. V. Kosynkin, A. L. Higginbotham, A. Sinitskii, J. R.
10 S. G. Sharapov, V. P. Gusynin, and H. Beck, Phys. Rev. Lomeda, A.Dimiev,B. K.Price, J.M. Tour,Nature458,
B 69, 075104 (2004). 872 (2009).
11 H. Fukuyama,J. Phys.Soc. Jpn.76, 043711 (2007). 41 L. Y. Jiao, L. Zhang, X. R. Wang, G. Diankov,H. J. Dai,
12 M. Nakamura,Phys. Rev.B 76, 113301 (2007). Nature458, 877 (2009).
13 M.KoshinoandT.Ando,Phys.Rev.B75,235333(2007). 42 S. Iijima, Nature(London) 354, 56 (1991).
14 A.Ghosal,P.Goswami, andS.Chakravarty,Phys.Rev.B 43 S. Iijima, T. Ichihashi and Y. Ando, Nature 356, 776
75, 115123 (2007). (1992).
15 T. Ando,Physica E 40, 213 (2007). 44 J.Liu,Z.Ma,A.R.Wright,andandChaoZhang,J.Appl.
16 M. Koshino, Y. Arimura, and T. Ando, Phys. Rev. Lett. Phys. 103, 103711 (2008).
102, 177203 (2009). 45 C. Berger, Z. Song, X. Li, X. Wu,N. Brown, C. Naud,D.
17 M.KoshinoandT.Ando,Phys.Rev.B81,195431(2010). Mayou, T. Li, J. Hass, A. N. Marchenkov, E. H. Conrad,
18 S. A.Safran, Phys.Rev. B 30, 421 (1984). P.N.First,and W.A.deHeer,Science312, 1191 (2006).
19 M.KoshinoandT.Ando,Phys.Rev.B76,085425(2007). 46 J. Hass, R. Feng, J. E. Millan-Otoya, X. Li, M. Sprinkle,
20 S. A. Safran and F. J. DiSalvo, Phys. Rev. B 20, 4889 P. N. First, W. A. de Heer, E. H. Conrad, and C. Berger,
(1979). Phys. Rev.B 75, 214109 (2007).
21 J. Blinowski and C. Rigaux, J. Phys. (Paris) 45, 545 47 G. Li, A. Luican and E. Y. Andrei, Phys. Rev.Lett. 102,
(1984). 176804 (2009).
22 R.SaitoandH.Kamimura,Phys.Rev.B33,7218(1986). 48 J. M. B. Lopes dos Santos, N. M. R. Peres, A.H. Castro
23 H.AjikiandT. Ando,J.Phys.Soc.Jpn.62, 2470 (1993); Neto, Phys.Rev. Lett.99, 256802 (2007).
J. Phys.Soc. Jpn. 63, 4267 (1994) (Erratum). 49 S. Latil, V. Meunier, and L. Henrard, Phys. Rev. B 76,
24 H.AjikiandT. Ando,J.Phys.Soc.Jpn.64, 4382 (1995). 201402(R) (2007).
25 M. Yamamoto, M. Koshino, and T. Ando, J. Phys. Soc. 50 S. Shallcross, S. Sharma and O.A. Pankratov, Phys. Rev.
Jpn. 77, 084705 (2008). Lett. 101, 056803 (2008).
26 M. Nakamura and L. Hirasawa, Phys. Rev. B 77, 045429 51 E. J. Mele, Phys.Rev.B 81, 161405 (2010).
(2008). 52 N. H. Shon and T. Ando, J. Phys. Soc. Jpn. 67, 2421
27 A. H. Castro Neto, F. Guinea, N. M. Peres, K. S. (1998).
Novoselov, and A. K. Geim, Rev. Mod. Phys. 81, 109 53 T. Ando, J. Phys. Soc. Jpn. 74, 777 (2005).
(2009). 54 V. P. Gusynin and S. G. Sharapov, Phys. Rev. Lett. 95,
28 A. Kobayashi, Y. Suzumura, and H. Fukuyama, J. Phys. 146801 (2005).
Soc. Jpn. 77, 064718 (2008). Hall effect and orbital dia- 55 N. M. R. Peres, F. Guinea, and A. H. Castro Neto, Phys.
magnetism in zerogap state of molecular conductor α- Rev.B 73, 125411 (2006).
(BEDT-TTF)2I3 56 P. Marconcini, and M. Macucci, Rivista Del Nuovo Ci-
29 M.Fujita,K.Wakabayashi,K.Nakada,andK.Kusakabe, mento 34 489 (2011).
J. Phys.Soc. Jpn. 65, 1920 (1996). 57 K. S,Krishnan and N.Ganguli, Nature 139, 155 (1937).
30 K. Nakada, M. Fujita, G. Dresselhaus, and M. S. Dressel-