Many-body effects on graphene conductivity: Quantum Monte Carlo calculations

0.36 MB·

by  D. L. Boyda

Many-body effects on graphene conductivity: Quantum Monte Carlo calculations
D. L. Boyda,1,2,∗ V. V. Braguta,2,3,4,† M. I. Katsnelson,5,6,‡ and M. V. Ulybyshev7,8,§
1Far Eastern Federal University, Vladivostok, 690091 Russia
2Institute of Theoretical and Experimental Physics, 117259 Moscow, Russia
3School of Biomedicine, Far Eastern Federal University, Sukhanova 8, Vladivostok, 690950 Russia
4Moscow Institute of Physics and Technology, Institutskii per. 9, Dolgoprudny, Moscow Region, 141700 Russia
5Radboud University, Institute for Molecules and Materials,
Heyendaalseweg 135, NL-6525AJ Nijmegen, The Netherlands
6UralFederalUniversity,TheoreticalPhysicsandAppliedMathematicsDepartment, MiraStr. 19,620002Ekaterinburg, Russia
6 7Institute of Theoretical Physics, University of Regensburg,
1 D-93053 Germany, Regensburg, Universitatsstrasse 31
0 8Institute for Theoretical Problems of Microphysics,
2 Moscow State University, Moscow, 119899 Russia

g
Optical conductivity of graphene is studied using Quantum Monte Carlo calculations. We start
u
fromEuclideancurrent-currentcorrelatorandextractσ(ω)fromGreen-KuborelationsusingBackus-
A
Gilbert method. Calculations were performed both for long-range interactions and taking into
account only contact term. In both cases we vary interaction strength and study its influence on
2
optical conductivity. Wecompare ourresults with previoustheoretical calculations choosing ω≈κ

] thusworkingintheregionoftheplateauinσ(ω)whichcorrespondstoopticalconductivityofDirac
l quasiparticles. No dependence of optical conductivity on interaction strength is observed unless
e
we approach antiferromagnetic phase transition in case of artificially enhanced contact term. Our
-
r results strongly support previous theoretical studies claimed very weak regularization of graphene
t
s conductivity.
.
t
a PACSnumbers: 73.22.Pr,05.10.Ln,11.15.Ha
m Keywords: grapheneconductivity, Coulombinteraction, Monte-Carlocalculations
-
d
n 1. INTRODUCTION onto p subspace results in a noticeable screening of ef-
z
o fective interaction potential at intermediate distances in
c Low-energyelectronicpropertiesofgrapheneinsingle- comparison with the bare Coulomb [10] which does not
[
particle approximation are determined by fermionic ex- change qualitatively perturbative results [11] but shifts

2 citations in p band with a linear relativistic dispersion thepointofsemimetal-insulatortransitionmakingfreely
z
v relation = v p and the Fermi velocity v c/300 suspended graphene semimetallic [12].
5 E F| | F ≃
1 [1–5]. Due to the smallness of the Fermi velocity vF ≪c An interesting observable for which one can expect
3 magneticinteractionbetweenfermionexcitationsandre- sizable renormalization is the optical conductivity σ(ω).
5 tardation effects are negligible. As a result, the interac- Theσ(ω) describeselectriccurrentingrapheneresulting
0 tionbetweenquasiparticlesisinstantaneousCoulomblaw fromexternalelectricfieldj(ω)=σ(ω)E(ω). Inthenon-
1. with effective coupling constant interacting theory of Dirac quasiparticles the conductiv-
0 ity does not depend on frequency in the limit of zero
1 c 1
6 α= (1) temperature and zero doping. It equals to[14] σ = 1
1 ǫvF 137 in the units of e2. The difference of the σ(ω) fro0m th4e
: ~
v where ǫ is dielectric permittivity of surrounding media. value σ can be attributed to the interaction between
0
i
X For suspended graphene the effective coupling constant quasiparticles. Experimental data on the optical con-
is sufficiently large α 2. Quite strong interaction be- ductivity did not find [15–17] any deviation from the re-
r 0
≈
a tween fermion excitations in graphene results in a rich sult of single-particle theory, within experimental error
variety ofphenomena [6]. In particular,“relativistic”in- bar, a few percent. This seems to be surprising keep-
variance of the theory is destroyed by renormalization ing in mind that many-body effects change dramatically
of the Fermi velocity which turns out to be strongly thesingle-electronspectrumandelectroncompressibility
k-dependent, k is the wave vector. This prediction [7] [8, 9].
has been recently confirmedexperimentally [8, 9]. Accu-
There are a lot of papers devoted to the leading-order
rate procedure of mapping of full electron Hamiltonian
perturbative correction to the value σ [18–26] which
0
present different results. In general, all these papers can
be divided into two groups. First of all, it was provenin
Ref. [23]thatnon-renormalizationofconductivity inthe
∗Electronicaddress: boyda˙[email protected]
limit of zero frequencies for weak on-site interaction is a
†Electronicaddress: [email protected]
‡Electronicaddress: [email protected] rigorous result following from symmetry required exact
§Electronicaddress: [email protected] cancellationofself-energyandvertexcontributions. This
2
work was based on renormalization group treatment of 2. MODEL DESCRIPTION AND
Hubbard model on honeycomb lattice. Thus, long-range CALCULATION OF CURRENT-CURRENT
interactionswerenottakenintoaccount. Earlier,theab- CORRELATOR
sence of renormalization was demonstrated phenomeno-
logically based on Fermi-liquid theory [24]. Whereas the Our model Hamiltonian consists of tight-binding term
system of electrons with contact interaction at the hon- andinteractionpartwhichdescribesthefullelectrostatic
eycomb lattice is probably Fermi liguid, for the case of interaction between quasiparticles :
long-rangeinteractionsthisisnotthecase[6]. Inthesec-
otanidl ognrouopptoicfaplacpoenrsduthcetivinitflyueonfcmeoafsslloensgs-rDainrgaecCfeorumloiomnbs Hˆ =−κ<Xx,y>(cid:16)aˆ†y,↑aˆx,↑+aˆ†y,↓aˆx,↓+h.c.(cid:17)
wasstudied. Therenormalizationofopticalconductivity + m(aˆ† aˆ aˆ† aˆ )
was found: x,↑ x,↑− x,↓ x,↓
X
x={1,ξ}
m(aˆ† aˆ aˆ† aˆ )
σσ(ω) =1+Cαeff +O(α2eff), (2) −x=X{2,ξ} x,↑ x,↑− x,↓ x,↓
0
1
+ V qˆ qˆ , (3)
xy x y
2
with αeff logarithmically dependent on frequency, but Xx,y
thevalueoftheconstantC wasdisputable. Sincealmost
whereκ=2.7eVishoppingbetweennearest-neighbours,
allofthesepaperswerebasedonDiracconesapproxima- aˆ† , aˆ and aˆ† , aˆ are creation/annihilation opera-
tion, C was essentially dependent on the regularization x,↑ x,↑ x,↓ x,↓
tors for spin up and spin down electrons at π-orbitals.
schemeusedincalculations. TheauthorsofRefs. [20,22]
Spatial index x = s,ξ consists of sublattice index
claimed rather large renormalization with C = 0.2...0.5. { }
s = 1,2 and two-dimensional coordinate ξ = ξ ,ξ of
On the other hand, Mishchenko [18] and Kotikov [19] { 1 2}
theunitcellinrhombiclattice. Periodicalboundarycon-
observed very small constant C 0.01. The authors of
≈ ditions areimposed in both spatial directions. The mass
Ref. [25] obtained C 0.26 working on honeycomb lat-
≈ term has different sign at different sublattices. Accord-
ticethususinganaturalcutoff. Inthemostrecentpaper
ing to our algorithm, we should add it to eliminate zero
[26] the authors calculated C taking into account non-
modes from fermionic determinant. Of course it means
zero size of atomic orbitals in interaction term. Similar
that we should check the dependence of final results on
to [18] they obtained very small renormalization.
the value of mass and this check is indeed performed in
This short review shows that renormalization of con- subsequent calculations.
ductivity in graphene is still an open, and quite contro- The matrix V is electrostatic interaction poten-
xy
versial, theoretical question. It seems that experiment tial between sites with coordinates x and y and qˆ =
x
supports small corrections to conductivity [15–17], thus aˆ† aˆ +aˆ† aˆ 1istheoperatorofelectricchargeat
there should exist a considerable cancellation between laxt,t↑icxe,↑site xx,.↓ Wx,↓e−use two sets of interaction potentials.
self-energy and vertex corrections. From the other side, The first one (“long range interaction”) consists of phe-
therearenotanyclearphysicalreasonswhichstaybehind nomenological potentials calculated by the constrained
thiscancellation. Inparticular,thisisnotthecaseforthe RPA method [10] at small distances. It continues at
other response function, electron compressibility, where large distances by ordinary Coulomb V 1/r . In
xy xy
self-energycorrectionsaredominant[6,9]. Whatismore general this set up corresponds to suspende∼d graphene,
important it is not clear if this cancellation takes place further details can be found in Ref. [12]. In the second
for higher order perturbation theory. So, the absence set up (“short range interaction”) we keep only on-site
or considerable suppression of the many-body effects in interactionequal to 9.3 eV. This number, again, accords
optical conductivity is important theoretical puzzle. toRef. [10]. Inordertostudydependenceoninteraction
In this paper we are going to address the question strength in both cases we used rescaling of potentials
of many-body effects in the optical conductivity using where all of them were uniformly divided by dielectric
Quantum Monte-Carlo simulation which fully accounts permittivity of surrounding media ǫ. Thus suspended
interactions between quasiparticles in a nonperturbative graphene corresponds to the point ǫ=1 on the plots for
way. In the second section we describe our model and long range interaction.
give a short review of Quantum Monte Carlo algorithm. All calculations were performed using Hybrid Monte-
Since the calculation of optical conductivity is based on Carlo algorithm. Details of the algorithm are described
linearresponsetheory,wedescribealsothecalculationof inthepapers[12,27,28]. ThemethodisbasedonSuzuki-
current-current correlator. The third section is devoted Trotter decomposition. Partition function exp( βHˆ) is
−
to the analytical continuation methods employed to ob- represented in the form of a functional integral in Eu-
tain real-time conductivity σ(ω). We present the final clidean time. Inverse temperature is equal to number
dependence of conductivity on interaction strength both of time slices multiplied by the step in Euclidean time:
for model with long-range and short-range interaction δτN =β =1/T. Sincethealgorithmneedsforfermionic
t
and give a short discussion of the results in conclusion. fieldstobeintegratedout,weeliminateallfour-fermionic
3
1.1 terms in full Hamiltonian using Hubbard-Stratonovich
m=0.0 eV
m=0.05 eV transformation. The final form of the partition function
1.05 m=0.07 eV can be written as:
m=0.1 eV
m=0.15 eV
τ) 1 m=0.2 eV
G (0 0.95 Tre−βHˆ = ϕ e−S[ϕx,n] detM[ϕ ]2. (4)
τ)/ ∼Z D x,n | x,n |
G(
0.9
0.85 ϕx,nistheHubbard-Stratonovichfieldfortimeslicenand
spatial coordinate x. Particular form of fermionic oper-
0.8
ator M is described in Ref. [12]. The absence of sign
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
τ T problem (appearance of the squared modulus of the de-
terminant) is guaranteed by the particle-hole symmetry
FIG. 1: The ratio of current-current correlator
ingrapheneatneutralitypoint. Actionforhubbardfield
calculated in presence of interaction and for free
S[ϕ ] is also positively defined quadratic form for all
x,n
fermions (G(τ)/G (τ)) as a function of Euclidean time
0 variants of electron-electron interaction used in our pa-
τ for different bare masses. Temperature is equal to
per. Thus we can generate configurations of ϕ by the
T =0.5 eV and lattice size is 242 20. Interaction x,n
Monte-Carlo method using the weight (4) and calculate
×
strength corresponds to suspended graphene.
physicalquantities as averagesovergeneratedconfigura-
tions.
ε = 9.00 We extract conductivity from Euclidean correlation
1.04 ε = 3.00 function of electromagnetic currents:
ε = 1.80
ε = 1.28
1.02 ε = 1.00
τG ()0 1 ε = 0.81 G(τ)= Tbc Tr(cid:16)e−βHˆJˆb(ξ)e−τHˆJˆc(ξ)eτHˆ(cid:17), (5)
τG()/ 3√3Ns Xξ Tre−βHˆ
0.98
0.96 whereNs is the number ofunit cells inoursample,sum-
mationover repeatedindexes is implied. Jˆ(ξ) is the op-
b
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
eratorofelectromagneticcurrentflowingfromthesite in
τ T
the first sublattice with coordinate x = 1,ξ towards
(a) The plot corresponds to themodel with long-range { }
one of its nearest neighbours. There are three near-
interaction.
est neighbours and 3 possible directions of the current:
1.1 b=1,2,3. Matrix Tbc is defined as follows:
ε = 9.00 ε = 1.80 ε = 1.00
ε = 3.00 ε = 1.28 ε = 0.81
1.05
1, b=c
τ) 1 Tbc =d2(cid:26) 1/2, b=c , (6)
(0 − 6
G 0.95
τ)/
G( 0.9 where d = 0.124nm is the distance between nearest
neighbours in graphene lattice. We need in calculation
0.85 the explicit form of electromagnetic current operator:
0.8
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
τ T Jˆ(ξ)= Jˆ (ξ)=iκ aˆ† aˆ +h.c., (7)
b σ,b σ,x σ,y
(b) Theplot corresponds to themodel with pureon-site σX=↑,↓ σX=↑,↓
interaction.
where x = 1,ξ , y = 2,ξ+ρ and shifts in rhombic
FIG. 2: The ratio of current-current correlator { } { b}
lattice are defined as follows: ρ = 0,0 , ρ = 1,1 ,
calculated in presence of interaction (G(τ)) and the one 1 { } 2 {− }
ρ = 1,0 . The final form of the current-current cor-
for free fermions (G0(τ)) as a function of Euclidean re3lato{r−used}in Monte Caro procedure can be written
time τ for different interaction strengths. Temperature
in terms of some combinations of the lattice fermionic
is equal to T =0.5 eV and lattice size is 242 20. Both
× propagator Γ(x,τ,x,τ′) = M−1(x,τ,x,τ′) (see [27] for
plots show the data for zero bare mass.
derivation):
4
2T
bc
G(τ)= Re j (y ,x )Γ(x ,0,x ,τ)j (x ,y )Γ(y ,τ,y ,0) +
b 1 1 1 2 c 2 2 2 1
−3√3N h i
s x1,xX2,y1,y2(cid:2) (cid:3)
4T
bc
Re j (y,x)Γ(x,0,y,0) Re j (y,x)Γ(x,τ,y,τ) . (8)
b c
3√3N h · i
s Xx,y(cid:2) (cid:3) Xx,y(cid:2) (cid:3)
0.4
Here j is current vertex:
b
0.35
Jˆ (ξ)= aˆ† aˆ j (x,y), (9)
σ,b σ,x σ,y b 0.3
Xξ x=X{s,ξ} )
y={s′,ξ′} 2-h/ 0.25
and ... is average over configurations of hubbard field (e 0.2
generhatied according to the statistical weight (4). -σω() 0.15
The firstterminexpression(8)is namedas connected
0.1
part and the second term is named as disconnected part
0.05
of the correlation function G(τ). Connected part was
calculated using standard stochastic estimator approach 0
0 1 2 3 4 5 6
[30]. The disconnected part is very noisy thus leading ω/κ
to very large uncertainty in the Monte-Carlosimulation.
(a) σ¯(ω) extracted from Green-Kubo relations. Vertical error
Fortunately one can argue that it is small as compared
baris statistical uncertainty. Horizontal error bar shows the
to the connected part. For τ = 0 the disconnected part
width of corresponding resolution function. It is calculated as
6
is the average of the operators which are located at dif- a width of thepeak at half of its maximum height. Central
ferent Euclidean time slices (see equation (8)). To get point is placed at ωc: maximum of corresponding resolution
nonzero result for the disconnected contribution the op- function.
erators at different τ–slices must interact. However, the
pdcinoiotsnwecnroeaenrccnotteefiocdtntheoebdneFedteswiraemagenrinadvmdeiwlffsoeecisricteyvan.entrTτys–hassufmlesilcyatehlslneiiesncgsolucneopctmrtpirbpteuhastersiiemosdonnbionyfwottihhntheee ω) 3456 ωωωωωω000000======011223......516273617283 κκκκκκω,,,,,, 0ωωωωωω=cccccc0======, 001223ω......495361c=928055 0κκκκκκ,
correlation function (8). Explicit check of this fact was ,0
ω
done in the paper [27]. In particular, we refer to the the δ( 2
figure 11. Profile of electron-electron interaction in our
1
paper differs from the one used in [27], so we can not
0
compare results with the same ε directly. Nevertheless,
it’spossibletousetheantiferromagneticphasetransition -1
0 1 2 3 4 5 6
as a reference point. Since our setup of electron-electron ω/κ
interactioncorrespondstothepointsinthesemi-metallic
phase away from the phase transition, we can use data (b) Resolution functions δ(ω0,ω).
sets also for semi-metallic phase (ε = 10 and ε = 4) at
FIG. 3: Results of test analytical continuation for free
the figure 11 in [27] as an illustration of the fact that
fermions: σ¯(ω) and corresponding resolution functions.
disconnected part of current-current correlator is much
Regularization is absent: λ=1. Bare mass term is
less than the connected one. Actually, it’s even impossi-
introduced: m=0.05 eV. κ=2.7 eV is hopping
bletodistinguishdisconnectedpartfromzerobecauseof
between nearest neighbours.
statistical errors.
Technically,wecarryoutthemeasurementsofG(τ) at
lattice with N = 242 cells and 20 Euclidean timeslices
s
(δτ = 0.1eV−1), which corresponds to the temperature and finite discretization errors are of order of statistical
T =0.5eV.Thedataareobtainedforthesetof5fermion uncertainty 0.3%.
∼
massesm=0.05,0.07,0.1,0.15,0.2eV.Inordertostudy The bare data for current-current correlator in sus-
systematic errors appeared due to the finite volume and pended graphene is presented in fig. 1. This is the ra-
finite discretization step in Euclidean time direction we tio G(τ)/G (τ), where G (τ) is the correlatorcalculated
0 0
also measured the correlation function G(τ) at lattice for free fermions in the massless limit. Data for zero
with N = 182, N = 20 and lattice with N = 242, bare mass is obtained via extrapolation by the function
s t s
N =40(δτ is two times smaller in the latter case). The G (τ) = A(τ) +B(τ) m2 + C(τ) m4. It should be
t m
· ·
results of the measurements show that the finite volume notedthatrenormalizationofthe G(τ) stronglydepends
5
0.4
via Green-Kubo relations:
0.35
∞
0.3 G(τ)= σ(ω)K(τ,ω)dω, (10)
Z
) 0
2-hω) (e/ 00.2.25 K(τ,ω)= ωcoπsshi(nωh((βω/β2/−2)τ)). (11)
-σ( 0.15
We used Backus-Gilbert method, recently adopted for
0.1
similar tasks in lattice quantum chromodynamics [29].
0.05
According to this method we obtain convolutionof opti-
0 cal conductivity with resolution functions δ(ω ,ω):
0 1 2 3 4 5 6 0
ω/κ ∞
(a) σ¯(ω) extracted from Green-Kuborelations. Vertical error σ¯(ω0)= δ(ω0,ω)σ(ω)dω. (12)
Z
baris statistical uncertainty. Horizontal error bar shows the 0
width of corresponding resolution function. It is calculated as Resolution functions are defined as linear combinations:
a width of thepeak at half of its maximum height. Central
point is placed at ωc: maximum of corresponding resolution Nt/2
function. δ(ω ,ω)= q (ω )K(jδτ,ω), (13)
0 j 0
56 ωωω000===123...123123 κκκω,,, 0ωωω=ccc0===, 123ω...013c=636 0κκκ, wimhiezraeticooneffiofcitehnetswqijd(tωhXj0=)o1fcatnhebereosobltuatiinoendffurnomctiothne. mTihne-
4
widthisdefinedasD = ∞(ω ω )2δ(ω ,ω)dω. Thisre-
ω) 3 0 − 0 0
,0 lationguaranteesthe miRnimum widthof resolutionfunc-
ω
δ( 2 tion with center close (but not exactly equal to) ω0.
We also take into account the normalization condition
1
∞
δ(ω ,ω)=1. The final formula for conductivity and
0 0 0
cRoefficients in (13) has the following form (see [29] for
-1 details):
0 1 2 3 4 5 6
ω/κ
Nt/2
(b) Resolution functions δ(ω0,ω). σ¯(ω )= q (ω )G(jδτ), (14)
0 j 0
Xj=1
FIG. 4: Real calculation for suspended graphene
W(ω )−1R
(long-range interaction with ǫ=1). In this case 0 jk k
q (ω )= , (15)
regularizationis imposed: λ=1 4 10−5. Bare mass j 0 R W(ω )−1R
n 0 nm m
− ×
term is introduced: m=0.05 eV. κ=2.7 eV is hopping
between nearest neighbours. whereW(ω0)jk =λ 0∞∞(ω−ω0)2K(jδτ,ω)K(kδτ,ω)dω+
(1 λ)S and R =R K(nδτ,ω)dω. S is covariance
− jk n 0 jk
matrix of current-curRrent correlator. Summation over
repeated indices is implied. Parameter λ is used to reg-
on the fermion bare mass m. The smaller the fermion
ularize poorly-defined matrix W . It rules the width of
mass used in the calculation the smaller the change of jk
resolutionfunctions. Ifqualityofdataisgoodenough,we
the correlation function G(τ) due to the interaction.
can work without any regularization (λ= 1) and obtain
TheratioG(τ)/G (τ)fordifferentvaluesofinteraction
0 minimalwidthofthefunctionsδ(ω ,ω). Withdecreasing
0
isshowninfig. 2forthemasslesslimit. Onecanseethat
ofλ the width of resolutionfunctions increases but algo-
renormalization of the G(τ) is not very large. It doesn’t
rithmbecomesmoretoleranttostatisticalerrorsinG(τ).
exceed 5% for long range interaction. It becomes notice-
The examples of resolution functions are plotted in the
ableonlyforthe largestvaluesofshort-rangeinteraction
figure 3b for λ = 1 (without regularization) and in the
whichareclosetoantiferromagneticphasetransition(see figure4bfor1 λ=1 10−5. Valueofω givesthecoor-
c
below). − ×
dinateofmaximumofcorrespondingresolutionfunction.
One can see that it doesn’t coincide with ω precisely,
0
but is very close to it. ω is used as x-coordinate in the
c
plots for σ(ω).
3. ANALYTICAL CONTINUATION AND
Results of test calculation of conductivity for free
CALCULATION OF CONDUCTIVITY
fermions are shown in the figure 3a. Even in calculation
of current-current correlator for free fermions we used
ThenextstageistransformationofEuclideancurrent- thesamenumericalprocedure(stochasticestimator)[30]
current correlator to real-frequency optical conductivity as for real calculation in interacting model. Thus we
6
0.34 0.6 1.7 -σ(ωc), ωc=1.06 κ
) 0.5 1.6 prediction with C=0.26
-h 0.4 prediction with C=0.01
-2he/) 00.3.32 2σω-() (e/c-0000....123010 1 2 3 4 5 6 7 8 9 10 ωσ)/c0 111...345
) ( (1.0- )x105 -σ( 1.2
ωc 0.28
σ-( 1.1
1
0.26
ω=1.06κ
c 0.9
0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3
0.24 1/ε
0 0.05 0.1 0.15 0.2
m, eV FIG. 6: Dependence of σ¯(ω ) on interaction
c |ωc=1.06κ
strength in the case of long-range interactions. All
λF=IG1. 5:4Dep1e0n−d5e.nIcnesoetf:σ¯D(ωecp)e|nωcd=e1n.c06eκoofncobnadruecmtivaistsyfoorn calculations were performed for λ=1−4×10−5 and
m=0.05 eV. Two variants of analytical predictions for
− ×
regularizationconstant λ for fixed mass m=0.05 eV.
graphene conductivity renormalization(see eq. 2) are
Both plots are for long-range interactions with ǫ=1
plotted for reference.
(suspended graphene).
1.05 σ
σ ,(cid:0)(cid:1)(cid:2)(cid:3)(cid:4)(cid:5)(cid:6)(cid:4)(cid:1)(cid:5)(cid:6)(cid:7)(cid:8)(cid:9)(cid:5)(cid:4)(cid:0)(cid:1)
havesomeresidualstatisticaluncertaintieswhichgiveus ,(cid:10)(cid:0)(cid:1)(cid:11)(cid:2)(cid:7)(cid:8)(cid:1)(cid:11)(cid:6)(cid:4)(cid:1)(cid:5)(cid:6)(cid:7)(cid:8)(cid:9)(cid:5)(cid:4)(cid:0)(cid:1)
1
a possibility to check stability of the method against er-
rors in input data. It appears that test calculation for
0.95
0
free fermions can be performed even without any regu- σ
σ
larization, so λ=1 for all data in the figure 3. 0.9
Withthistestcalculationweshouldreproduceanalyti-
calresultforσ(ω)onhoneycomblatticeandthefigure3a 0.85
indeed demonstrates all essential properties of free spec-
0.8
tralfunctionσ(ω): transportpeakatzerofrequency,then 0 0.2 0.4 0.6 0.8 1 1.2 1.4
plateau around σ , then van Hove singularity and final ε
0
decay at large frequencies. We conclude that we should
FIG. 7: Left scale: dependence of conductivity σ˜
take resolution function with maximum around ω κ
c ≈ (calculated from the middle point of current-current
to reproduce right value of optical conductivity. This ω
c correlator,see eq. 16) on the interaction strength both
is large enough to get rid of the gap induced by finite
for long-range and short-range interactions. Phase
mass (we do not see any dependence on bare mass, see
transition into antiferromagnetic state is situated at
below) and small enough for results not to be influenced
ε=0.6...0.7 in the case of long-range interactions (see
by van Hove singularity.
[12]) and at ε=0.91 (see [13]) in the case of pure
In real Monte-Carlo calculations statistical errors are
on-site interaction.
sufficientlylargerandweneedsomeregularizationλ=1.
6
Figure4showsopticalconductivityσ¯(ω)andcorrespond-
ing resolution functions for suspended graphene (ǫ = 1,
long-range interaction). Since we are forced to intro- Again the dependence on λ is flat unless the instability
duce regularization in this case, the width of δ(ω ,ω) of analytical continuation causes sharp growth of errors
0
functions for large ω increases and we can not catch at λ 1. According to results of this test we choose
0 →
van Hove singularity in our data. Nevertheless, width of regularizationconstant λ=1 4 10−5.
− ×
δ(ω ,ω) function with maximum at ω κ is rather sta- Using this set up we plot the dependence of conduc-
0 c
≈
ble and we can still use it to estimate the value of σ(ω) tivity on interaction strength in the case of long-range
at plateau which corresponds to optical conductivity of interaction(fig. 6). Forreferenceweplotatthesamefig-
Diracquasiparticles. Thisresolutionfunctionwithmaxi- ure also two variants of analytical predictions according
mumatω =1.06κisusedinallsubsequentcalculations. to eq. 2 for two different scenarios: large renormaliza-
c
Now we should check independence of results on bare tion with C =0.26 and small renormalization C =0.01.
massintheHamiltonian(3)andregularizationλ. Figure Our data is consistent only with small renormalization
5demonstratesresultsofthischeckinthecaseoffullcal- scenariowhilelargerenormalizationofgrapheneconduc-
culation in interacting model. Main plot shows that we tivity is completely ruled out by our calculations.
don’t have any definite dependence on bare mass. Thus In order to estimate conductivity in the case of short-
we simply use the lowestpossible m=0.05eV in allcal- rangeinteractions,we used alternativeapproach. Unfor-
culations. Effect of regularization is shown in the inset. tunately, statistical errors in current-current correlator
7
are too large in this case and we are unable to extract to proximity to the AFM phase transition. Thus, our
σ(ω) directly from Green-Kubo relations. We used es- calculations are accurate enough to agree with analyti-
timate of optical conductivity from the middle point of calresults[23];beyondthephasetransition,perturbative
correlator (taken in the massless limit): arguments of that work are no more valid.
β2
σ˜ = G(β/2) . (16)
m→0
π | 4. CONCLUSIONS
The method was proposed in [27]. Physically it corre-
To conclude, straightforward Quantum Monte Carlo
sponds to conductivity smeared over some region of fre-
calculations of optical conductivity for interacting elec-
quencies near ω = 0. Width of the region is approxi-
tronsonhoneycomblatticeshowaveryweakdependence
mately equal to 2T:
on the coupling constant. Numerical data completely
β2 ∞ dω 2ω rulesoutthe“strongrenormalization”scenario. Itcanbe
σ˜ = π Z 2π sinh 1βω σ(ω). (17) consistentwithboth anaccidental“almostcancellation”
0 2 of self-energy and vertex corrections leading to a small
(cid:0) (cid:1)
coefficient C in Eq.(2) and an exact calculation due to
In the case of free fermions this approach gives us the
value σ˜ 0.28 which is also rather close to σ . Re- Wardidentity,similartothecaseofHubbardmodel[23].
0
≈ It is not clear yet, however, why such exact cancellation
sults for conductivity calculation from the middle point
ofcurrent-currentcorrelator(σ˜)areshowninthefigure7 should take place for long-range interactions. We hope
that computational results presented here will stimulate
in the case of short- and long-range interactions. Again,
this method doesn’t show any sufficient renormalization further analytical studies of this very complicated and
controversialissue.
of optical conductivity in the case of long-range inter-
action. In the case of short-range interactions we see Acknowledgements
decrease of conductivity for large interaction strength.
Stronger renormalization in the case of pure short The work of MU was supported by the DFG Grant
range interaction is connected to the fact that antifer- BU 2626/2-1 and by Grant RFBR-14-02-01261-a. MIK
romagnetic(AFM)phasetransitionisclosertotheε=1 acknowledges financial support from NWO via Spinoza
point in the case of short-range interaction. While it Prize. The work of DLB was supported by RFBR grant
takes place at ε = 0.6...0.7 in the model with long- 16-32-00362-mol-a. The work of VVB, which consisted
range Coulomb tail [12], pure on-site interaction causes of the calculation of the current-current correlator, was
phase transition at V = 3.78κ [13] which corresponds supportedbygrantfromtheRussianScienceFoundation
00
to ε=0.91inournotation(see the descriptionofthe in- (project number 16-12-10059). The authors are grateful
teractionsaftereq. (3)). Onecanconcludethatdecrease to FAIR-ITEP supercomputer center and to supercom-
of σ˜ in the case of short-range interactions is connected puter center of Moscow State University.
[1] P.R. Wallace, Phys.Rev.71, 622 (1947). stein, M. I. Katsnelson, and S. Blugel, Phys. Rev. Lett.
[2] J. W. McClure, Phys. Rev.104, 666 (1956). 106, 236805 (2011).
[3] G. W. Semenoff, Phys. Rev.Lett. 53, 2449 (1984). [11] N.Yu.Astrakhantsev,V.V.Braguta,andM.I.Katsnel-
[4] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, son, Phys.Rev.B 92, 245105 (2015).
M. I. Katsnelson, I. V. Grigorieva, S. V. Dubonos, and [12] M. V. Ulybyshev, P. V. Buividovich, M. I. Katsnelson,
A.A. Firsov, Nature 438, 197 (2005). and M. I. Polikarpov, Phys. Rev. Lett. 111, 056801
[5] Y.Zhang,Y.-W.Tan,H.L.Stormer,andP.Kim,Nature (2013).
438, 201 (2005). [13] F. F. Assaad and I. F. Herbut, Phys. Rev. X. 3, 031010
[6] V. N. Kotov, B. Uchoa, V. M. Pereira, F. Guinea, and (2013).
A.H. Castro Neto, Rev.Mod. Phys. 84, 1067 (2012). [14] M.I.Katsnelson,Graphene: Carbon inTwoDimensions
[7] J.Gonzalez,F.Guinea,andM.A.H.Vozmediano,Nucl. (Cambridge University Press, 2012).
Phys.B 424, 595 (1994). [15] Z.Q.Li, E.A.Henriksen,Z.Jiang, Z.Hao, M.C. M.P.
[8] D.C. Elias, R. V.Gorbachev, A.S.Mayorov, S. V.Mo- Kim, H. L. Stormer, and D. N. Basov, Nature Phys. 4,
rozov,A.A.Zhukov,P.Blake,L.A.Ponomarenko,I.V. 532 (2008).
Grigorieva,K.S.Novoselov,F.Guinea,andA.K.Geim, [16] K.F.Mak,M.Y.Sfeir,Y.Wu,C.H.Lui,J.A.Misewich,
NaturePhys. 7, 701 (2011). and T. F. Heinz, Phys. Rev.Lett. 101, 196405 (2008).
[9] G. L. Yu,R. Jalil, B. Belle, A. S. Mayorov, P. Blake, F. [17] R.R.Nair,P.Blake,A.N.Grigorenko, K.S.Novoselov,
Schedin, S. V. Morozov, L. A. Ponomarenko, F. Chiap- T.J.Booth,T.Stauber,N.M.R.Peres,andA.K.Geim,
pini, S. Wiedmann, U. Zeitler, M. I. Katsnelson, A. K. Science 320, 1308 (2008).
Geim,K.S.Novoselov,andD.C.Elias,Proc.Nat.Acad. [18] E. G. Mishchenko, Europhys.Lett. 83, 17005 (2008).
Sci. 110, 3282 (2013). [19] S. Teber, A.V. Kotikov, Europhys. Lett. 107, 57001
[10] T. O. Wehling, E. Sasioglu, C. Friedrich, A. I. Lichten- (2014).
8
[20] I.F. Herbut, V. Juricic, and O. Vafek, Phys. Rev. Lett. [26] J.M. Link,P. P. Orth,D. E. SheehyandJ. Schmalian,
100, 046403 (2008). ArXiv:1511.05984.
[21] D.E.SheehyandJ.Schmalian,Phys.Rev.B80,193411 [27] P.V.BuividovichandM.I.Polikarpov,Phys.Rev.B86,
(2009). 245117 (2012).
[22] V. Juricic, O. Vafek, and I.F. Herbut, Phys. Rev. B 82, [28] D. Smith and L. von Smekal, Phys. Rev. B 89, 195429
235402 (2010). (2014).
[23] A. Giuliani, V. Mastropietro, and M. Porta, Phys. Rev. [29] B.B. Brandt,A. Francis,H. B.Meyer,andD. Robaina,
B 83, 195401 (2011). Phys. Rev.D 92, 094510 (2015).
[24] M. I. Katsnelson, Europhys. Lett.84, 37001 (2008). [30] T. Degrand, C. DeTar, Lattice Methods for Quantum
[25] B.Rosenstein,M.Lewkowicz,andT.Maniv,Phys.Rev. Chromodynamics (World Scientific Publishing, 2006).
Lett.110, 066602 (2013).