Velocity renormalization in graphene from lattice Monte Carlo

arXiv:1304.1711v2 [cond-mat.str-el] 26 Mar 2014
Velocity renormalization in graphene from lattice
Monte Carlo
Joaquín E. Druta and Timo A. Lähde∗b
a
Department of Physics and Astronomy, University of North Carolina, Chapel Hill, North
Carolina 27599-3255, USA
b Institute for Advanced Simulation, Institut für Kernphysik, and
Jülich Center for Hadron Physics, Forschungszentrum Jülich, D–52425 Jülich, Germany
E-mail: [email protected], [email protected]
We compute the Fermi velocity of the Dirac quasiparticles in clean graphene at the charge neutrality point for strong Coulomb coupling αg . We perform a Lattice Monte Carlo calculation
within the low-energy Dirac theory, which includes an instantaneous, long-range Coulomb interaction. We find a renormalized Fermi velocity vFR > vF , where vF ≃ c/300. Our results are
consistent with a momentum-independent vFR which increases approximately linearly with αg ,
although a logarithmic running with momentum cannot be excluded at present. At the predicted
critical coupling αgc for the semimetal-insulator transition due to excitonic pair formation, we
find vFR /vF ≃ 3.3, which we discuss in light of experimental findings for vFR /vF at the charge
neutrality point in ultra-clean suspended graphene.
PACS: 73.22.Pr, 73.63.Bd, 71.30.+h, 05.10.Ln
31st International Symposium on Lattice Field Theory LATTICE 2013
July 29 – August 3, 2013
Mainz, Germany
∗ Speaker.
c Copyright owned by the author(s) under the terms of the Creative Commons Attribution-NonCommercial-ShareAlike Licence.
http://pos.sissa.it/
Velocity renormalization in graphene from lattice Monte Carlo
Timo A. Lähde
1. Introduction
Graphene, a two-dimensional membrane of carbon atoms with unique electronic properties, is
characterized by a low-energy spectrum of Dirac quasiparticles, with a Fermi velocity vF ≃ 1/300
of the speed of light in vacuum [1, 2]. As the strength of the Coulomb interaction between
the quasiparticles is controlled by αg ≡ e2 /(4πε vF ) ≃ 2.2/ε , the role of interactions can be enhanced to the point that graphene may resemble Quantum Electrodynamics (QED) in a strongly
coupled regime [3]. In particular, the unscreened, long-range Coulomb interaction in graphene
leads to non-trivial velocity renormalization effects. At weak coupling, a logarithmic running
vF (n)/vF (n0 ) = 1 + (αg /4) ln(n0 /n) with carrier density n is found [4], such that vFR is expected
to become large in the vicinity of the Dirac point (n = 0). On the experimental side, logarithmic
velocity renormalization is most prominent in ultra-clean suspended graphene [6] and on boron
nitride (BN) substrates [7], with vFR /vF ≃ 2 − 3 in the former case, where ε = 1.
Electron-electron interactions may also trigger a semimetal-insulator transition due to excitonic pairing of quasiparticles and holes at a critical coupling αgc . For graphene, Lattice Monte
Carlo (LMC) has been applied to the Dirac theory using a contact Thirring interaction [8] and a
long-range Coulomb interaction [9 – 11], and to the tight-binding Hamiltonian using interactions
of the Hubbard [12 – 14] and Coulomb [15 – 17] types. For the Dirac theory, αgc ≃ 1.11 ± 0.06
was found [9], to be compared with αgc ≃ 0.9 ± 0.2 in the tight-binding + Coulomb approach [16].
While such a transition has not yet been observed, the empirical vFR /vF ≃ 2 − 3 in suspended
graphene is indicative of interaction-induced spectral changes [3, 6].
2. Lattice Monte Carlo
The LMC treatment of graphene uses the linearized low-energy Hamiltonian [18, 19] with an
instantaneous Coulomb interaction, such that Aµ ≡ (A0 ,~0). This gives the Euclidean (continuum)
action
Nf Z
Z
1
SE ≡ 2 d 3 x dt (∂i A0 )2 + ∑ d 2 x dt ψ¯ a D[A0 ]ψa ,
(2.1)
2g
a=1
where g2 ≡ e2 /ε , with ε ≡ (1 + κ )/2 for a substrate with dielectric constant κ . Here, ψa is a
four-component Dirac field in 2 + 1 dimensions with ψ¯ ≡ ψ † γ0 , A0 is the gauge (photon) field in
3 + 1 dimensions, and N f = 2 for a graphene monolayer. The Dirac operator is
D[A0 ] = γ0 (∂0 + iA0 ) + vF
2
∑ γk ∂k + m0,
(2.2)
k=1
where the γµ satisfy the Euclidean Clifford algebra {γµ , γν } = 2δµν , and the bare fermion mass
m0 provides an infrared regulator for modes that would be massless when the U(2N f ) chiral symmetry of Eq. (2.1) is spontaneously broken. The lattice version of Eq. (2.1) is formulated in terms
of “staggered” fermions, i.e. one-component Grassmann variables χ , χ¯ , which partially retain the
U(2N f ) chiral symmetry of Eq. (2.1) at finite lattice spacing (for staggered fermions in odd dimensions, see Ref. [20]). The fermionic part of Eq. (2.1) is given for N f = 2 in terms of staggered
fermions by ∑n,m χ¯ n Kn,m [θ0 ] χm , where n ≡ (n0 , n1 , n2 ) = (t, x, y) and m denote lattice sites on
2
Velocity renormalization in graphene from lattice Monte Carlo
Timo A. Lähde
a 2 + 1 dimensional fermion “brane”. This square space-time lattice embedded in a cubic lattice
should not be understood as a tight-binding theory. The staggered Dirac operator is
Kn,m [θ0 ] =
1
λ
†
(δn+e0 ,m U0,n − δn−e0 ,m U0,m
) + ∑ ηni (δn+ei ,m − δn−ei ,m ) + m0 δn,m ,
2a
2a i
(2.3)
where ηn1 = (−1)n0 , ηn2 = (−1)n0 +n1 , with ei a unit vector in lattice direction i. The invariance of
Eq. (2.1) under spatially uniform, time-dependent gauge transformations is retained by the gauge
links U0 ≡ exp(iθ0 ). We perform LMC calculations for vF = 1 (thus g2 → g2 /vF ) and a/ax = 1,
where a ≡ at is the temporal lattice spacing, and thus λ = 1. At non-zero αg , we have λR ≡
vFR (a/ax )R from which vFR /vF can be obtained, once the asymmetry (a/ax )R is known.
We compute the renormalized λR and mR from the staggered fermion propagator C f (t, x, y) ≡
−1 i on an N 2 × N space-time lattice with N /4 integer. Here n is
hχ (t, x, y)χ¯ (t0 , x0 , y0 )i = hKn,n
x
t
x,t
0
0
an arbitrary point of reference, and the brackets denote an average over ensembles of gauge field
configurations, obtained as a function of β ≡ vF /g2 = 1/(4παg ) and m0 , with the Hybrid Monte
Carlo algorithm. From C f (t, x, y), we form C f t (t, p1 , p2 ) ≡ ∑x,y exp(−ip·x)C f (t, x, y), for momenta
pi = 2π n/Nx , with n = −Nx /4, . . . , Nx /4, t = 0, . . . , Nt − 1, and summation over even lattice sites,
defined by (−1)t+x+y = 1. The expression for the “temporal correlator” CRft with renormalized mR ,
λR and wave function renormalization ZR (for a = 1) is [21]
CRft (t, p1 , p2 ) = ZR Gt (p1 , p2 ),
(2.4)
for t = 0, 2, . . . , Nt − 2, and
CRft (t, p1 , p2 ) = −
ZR
exp(iB0 ) Gt+1 (p1 , p2 ) − exp(−iB0 ) Gt−1 (p1 , p2 ) ,
2mR
(2.5)
for t = 1, 3, . . . , Nt − 1, with anti-periodic boundary conditions. The function Gt (p1 , p2 ) is
Gt (p1 , p2 ) ≡
N
C2 (µt ) − B2 (B0 )
× A(B0 )C(µt ) cos(B0 t ⋆ ) sinh(µt t ⋆ ) + B(B0 )D(µt ) sin(B0t ⋆ ) cosh(µt t ⋆ )
+ iA(B0 )C(µt ) sin(B0t ) sinh(µt t ) − iB(B0 )D(µt ) cos(B0t ) cosh(µt t ) ,
⋆
⋆
⋆
⋆
(2.6)
where A(x) ≡ cos(xNt /2), B(x) ≡ sin(xNt /2), C(x) ≡ cosh(xNt /2), D(x) ≡ sinh(xNt /2), t ⋆ ≡
Nt /2 − t, N ≡ 2mR / sinh(2µt ), and sinh2 (µt ) ≡ m2R + λR2 sin2 (p1 ) + λR2 sin2 (p2 ). This expression
for CRft includes a constant “background field” B0 ≡ hθ0 i, as hθ0 i 6= 0 in a finite volume. B0 is
roughly bounded by ±π /Nt [21], and the imaginary part of CRft vanishes for B0 → 0.
We also define the fermion “spatial correlator” C f x (x, ω , p2 ) ≡ ∑t,y exp(−ip · x)C f (t, x, y), for
ω = 2π (n− 1/2)/Nt (due to the anti-periodic temporal boundary conditions), n = −Nt /4, . . . , Nt /4,
spatial slices x = 0, . . . , Nx − 1, and summation over even lattice sites. The expression for CRfx
can be inferred from CRft . The function Gx (ω , p2 ) is obtained from Eq. (2.6) by first exchanging
sin ↔ cos, followed by t → x, Nt → Nx , µt → µx and B0 → 0. In addition, we have mR → mR /λR
in Eq. (2.5), mR → mR /λR2 in the expression for N, and sinh2 (µx ) ≡ m2R /λR2 + sin2 (ω + B0 )/λR2 +
sin2 (p2 ). Unlike Eqs. (2.4) and (2.5), CRfx is real-valued, with periodic boundary conditions.
3
Velocity renormalization in graphene from lattice Monte Carlo
Nt = 28, Nx = 28, Nb = 8
Cft, m0 = 0.020
0.010
0.010
0.3
0.005
3.25
0.0025
0.0025
3
2.75
2.5
2.5
2.25
2
2
1.75
1.75
0.15
1.5
0.05
0.2
0.1
0.15
β
Nt = 32, Nx = 32, Nb = 32
4
Cft, m0 = 0.015
3.75
0.010
Cft, m0 = 0.010
3.5
0.0075
3.5
0.005
0.005
3.25
0.004
0.0025
0.003
λR
3
2.75
2.5
2.25
2.25
2
2
1.5
0.05
0.2
0.15
0.15
0.1
0.1
0.05
0.05
0.0025
2.75
2.5
1.75
0.2
0.008
3.25
3
0.25
mR
β
4
0.25
0.2
Nt = 32, Nx = 32, Nb = 12
3.75
0.3
2.75
2.25
0.1
1/β = 5.0
6.0
7.0
8.0
9.0
10.0
11.0
12.0
13.0
14.0
15.0
mR
λR
3
λR
0.015
3.5
0.005
3.25
Cfx, m0 = 0.020
3.75
0.015
3.5
λR
1/β = 5.0
6.0
7.0
8.0
9.0
10.0
11.0
12.0
13.0
14.0
15.0
4
3.75
1.5
0.05
0.35
Nt = 28, Nx = 28, Nb = 8
4
Timo A. Lähde
1.75
0.1
β
0.15
1.5
0.05
0.1
0.15
β
0
0
0
0.005
0.01
m0
(a) λR as a function of β ≡ 1/(4παg ) and m0 , as obtained
from C f t and C f x on 283 × 8 (upper panels), 323 × 12 and
324 (lower panels) lattices. Note that the dependence on m0
is negligible, and that the results obtained from C f t and C f x
agree closely.
0.015
0.02
0
0.005
0.01
0.015
0.02
m0
(b) Left panel: mR as a function of β −1 ≡ 4παg and
m0 , obtained from C f t (colored connected symbols) and
C f x (black unconnected symbols) on a 283 × 8 lattice.
Right panel: mR obtained from C f t on 283 × 8 (colored
connected symbols) and 323 × 12 lattices (black unconnected symbols).
Figure 1: LMC calculation of λR and mR (for a = 1). All lines are intended as a guide to the eye, and
statistical errors are comparable to the size of the symbols.
3. Results
In Fig. 1(a), we show λR as obtained from a chi-square fit of Eqs. (2.4) and (2.5) to LMC
data. While Eq. (2.1) is gauge invariant, C f t and C f x are not, and thus a gauge fixing condition is
imposed on each configuration in order to obtain stable results. For C f t , a correlated fit is performed
for all (t, p1 , p2 ), and in the case of C f x for all (x, ω , p2 ). The fitted parameters are mR , λR , B0 and
ZR . Our lattices have Nt = Nx , and the length of the “bulk” dimension (where only the photons
propagate) is denoted Nb . We use the notation Nx3 × Nb , and simply Nx4 when Nb = Nt = Nx . On
the 283 × 8 lattice, LMC data is available for (inverse) lattice couplings 5.0 ≤ β −1 ≤ 15.0, and
for bare quasiparticle masses 0.0025 ≤ m0 ≤ 0.020, with slightly more restrictive data sets on the
323 × 12 and 324 lattices. We find that λR increases monotonically as a function of αg from the
non-interacting value of unity, with no appreciable differences between λR as obtained from C f t
and C f x . We find the dependence on m0 to be almost negligible. Finite size effects for λR are small,
and the fitted values of B0 agree closely with hθ0 i.
In Fig. 1(b), we show the physical quasiparticle mass mR as a function of β −1 and m0 , with
emphasis on asymmetries between the temporal and spatial correlations, and on finite size effects.
4
Velocity renormalization in graphene from lattice Monte Carlo
Timo A. Lähde
As Eq. (2.1) treats space and time asymmetrically, the spatial and temporal correlation lengths ξ
may exhibit unequal critical scaling, such that ξs ∝ |β − βc |−νs and ξt ∝ |β − βc |−νt . The dynamical
critical exponent z ≡ νt /νs is an important characteristic of a quantum critical point (QCP), and
implies that the dispersion relation is modified to E ∝ pz . At large N f , Ref. [5] predicted z ≃ 0.8
for graphene in the strong-coupling limit. However, arguments have also been put forward which
indicate z = 1 for a QCP with d < 4 in theories with a long-range Coulomb interaction [23]. From
Fig. 1(b), we find that the values of mR obtained from C f t and C f x agree very closely for β −1 ≤ 11.0,
which is consistent with z ≃ 1. We also find no sign of non-linear dispersion. A more accurate
analysis is possible following Ref. [8], in terms of the equation of state (EOS)
(δ ·βm −1)/νt
m0 = A(β − βc ) mR
δ ·βm /νt
+ B mR
,
(3.1)
for mR computed from C f t , where δ and βm are critical exponents characterizing the QCP at β = βc .
An analogous EOS with νt → νs can be given for mR as obtained from C f x . We also find from
Fig. 1(b) that finite size effects are under control for β −1 ≤ 9.0, but not for smaller β (stronger
coupling), especially in the region of the phase diagram where a dynamically generated gap exists,
and limm0 →0 mR 6= 0. In principle, Eq. (3.1) can be used to compute z and mR in the limit m0 → 0,
i.e. the gap in the quasiparticle spectrum. The lattice spacing asymmetry (a/ax )R could then be
inferred by measuring the gap in terms of temporal and spatial correlations. For reliable conclusions, such an EOS analysis should be combined with an extrapolation to infinite volume, similar
to that of Ref. [22] for QED4 . For the present analysis, we note from Fig. 1(b) that the values of mR
computed from C f t and C f x differ by no more than ≃ 10% at the smallest values of β . Assuming
this trend persists in the limits of infinite volume and vanishing m0 , we find 1.0 ≤ (a/ax )R ≤ 1.1
over the range of αg studied. In the absence of substantial indications for (a/ax )R 6= 1, we take λR
as a measure of vFR /vF .
4. Conclusions
In Fig. 2, we summarize our LMC results for vFR /vF as a function of αg ≡ 1/(4πβ ), and compare with available experimental data. Throughout our analysis, we have assumed that vFR is constant, while in reality we should expect it to run logarithmically with the momentum p and diverge
at the Dirac point, according to vF (p)/vF (p0 ) = 1 + (αg /4) ln(p0 /p), where p0 is the momentum
scale at which vFR = vF . At present, we cannot distinguish between this and the simpler expression
vFR /vF = 1 + Cαg , and much larger lattices appear to be required to detect a logarithmic running
of the Fermi velocity with p. On the other hand, we find a pronounced dependence of vFR /vF on
αg . We find that vFR increases linearly with αg from the non-interacting value vF up to αg ≃ 0.5,
above which the increase becomes more rapid. At the predicted critical coupling αgc ≃ 1.1, we
find vFR (αgc )/vF ≃ 3.3 within our present linearized treatment of the velocity renormalization.
Since all of the empirical vFR (p = 0)/vF fall short of this, we find a plausible explanation for the
non-observation of excitonic insulating phases in graphene monolayers. We note that the result of
Ref. [6] is tantalizingly close to ≃ 3.3, which suggests that further refinements in the quality of
suspended graphene may suffice to trigger the excitonic instability. It would be of interest to study
the logarithmic running of vFR with momentum on larger lattices.
5
Velocity renormalization in graphene from lattice Monte Carlo
Timo A. Lähde
4
suspended graphene
3.75
on BN substrate
3.5
3
LMC, 28 x 8
3.25
3
LMC, 32 x 12
vFR / vF
3
2.75
2.5
2.25
2
1.75
1.5
1.25
1
0.75
0
0.2
0.4
0.6
0.8
αg
1
1.2
1.4
Figure 2: Plot of vFR (αg )/vF for asymmetry (a/ax )R ≃ 1. The black line is vFR /vF = 1 + Cαg , and the
predicted insulating phase occurs for αg > αgc (gray shaded area). We find a linear dependence of vFR /vF
on αg up to αg ≃ 0.5 (solid black line). Note that vFR (αg = 0)/vF ≡ 1. Horizontal bands indicate empirical
data on vFR (p = 0)/vF from Ref. [6] (suspended graphene) and Ref. [7] (on BN substrate). Note that the
expected logarithmic momentum-dependence of vFR /vF cannot be resolved on present lattices.
Acknowledgments
We thank Lars Fritz, Simon Hands and Gerrit Schierholz for instructive discussions, and Jan
Bsaisou, Dean Lee and Ulf-G. Meißner for a careful reading of the manuscript. T. L. acknowledges
financial support from the Helmholtz Association (contract VH-VI-417). This work was supported
in part by an allocation of computing time from the Ohio Supercomputer Center (OSC).
References
[1] K. S. Novoselov et al., Electric Field Effect in Atomically Thin Carbon Films, Science 306, 666
(2004); Two-dimensional atomic crystals, Proc. Natl. Acad. Sci. U.S.A. 102, 10451 (2005);
Two-dimensional gas of massless Dirac fermions in graphene, Nature (London) 438, 197 (2005);
A. K. Geim and K. S. Novoselov, The rise of graphene, Nature Mater. 6, 183 (2007).
[2] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, The electronic
properties of graphene, Rev. Mod. Phys. 81, 109 (2009).
[3] V. N. Kotov, B. Uchoa, V. M. Pereira, F. Guinea, and A. H. Castro Neto, Electron-Electron
Interactions in Graphene: Current Status and Perspectives, Rev. Mod. Phys. 84, 1067 (2012).
[4] J. González, F. Guinea, and M. A. H. Vozmediano, Non-Fermi liquid behavior of electrons in the
half-filled honeycomb lattice (A renormalization group approach), Nucl. Phys. B424, 595 (1994).
[5] D. T. Son, Quantum critical point in graphene approached in the limit of infinitely strong Coulomb
interaction, Phys. Rev. B 75, 235423 (2007);
[6] D. C. Elias et al., Dirac cones reshaped by interaction effects in suspended graphene, Nature Phys. 7,
701 (2011).
[7] G. L. Yu et al., Interaction phenomena in graphene seen through quantum capacitance, Proc. Natl.
Acad. Sci. U.S.A. 110, 3285 (2013).
6
Velocity renormalization in graphene from lattice Monte Carlo
Timo A. Lähde
[8] S. J. Hands and C. G. Strouthos, Quantum critical behavior in a graphenelike model, Phys. Rev. B 78,
165423 (2008); W. Armour, S. Hands, and C. Strouthos, Monte Carlo simulation of the
semimetal-insulator phase transition in monolayer graphene, ibid. 81, 125105 (2010); Monte Carlo
simulation of monolayer graphene at nonzero temperature, ibid. 84, 075123 (2011).
[9] J. E. Drut and T. A. Lähde, Is Graphene in Vacuum an Insulator?, Phys. Rev. Lett. 102, 026802
(2009); Lattice field theory simulations of graphene, Phys. Rev. B 79, 241405 (2009); Critical
exponents of the semimetal-insulator transition in graphene: A Monte Carlo study, ibid. 79,
165425(R) (2009).
[10] P. V. Buividovich, E. V. Luschevskaya, O. V. Pavlovsky, M. I. Polikarpov, and M. V. Ulybyshev,
Numerical study of the conductivity of graphene monolayer within the effective field theory approach,
Phys. Rev. B 86, 045107 (2012).
[11] E. Shintani and T. Onogi, Chiral symmetry breaking in lattice QED model with fermion brane,
[arXiv:1203.1091 [hep-lat]].
[12] S. Sorella and E. Tosatti, Semi-Metal-Insulator Transition of the Hubbard Model in the Honeycomb
Lattice, Europhys. Lett. 19, 699 (1992).
[13] T. Paiva, R. T. Scalettar, W. Zheng, R. R. P. Singh, and J. Oitmaa, Ground-state and finite-temperature
signatures of quantum phase transitions in the half-filled Hubbard model on a honeycomb lattice,
Phys. Rev. B 72, 085123 (2005).
[14] Z. Y. Meng, T. C. Lang, S. Wessel, F. F. Assaad, and A. Muramatsu, Quantum spin liquid emerging in
two-dimensional correlated Dirac fermions, Nature (London) 464, 847 (2010).
[15] R. Brower, C. Rebbi, and D. Schaich, Hybrid Monte Carlo simulation on the graphene hexagonal
lattice, PoS(Lattice 2011), 056 (2011).
[16] P. V. Buividovich and M. I. Polikarpov, Monte Carlo study of the electron transport properties of
monolayer graphene within the tight-binding model, Phys. Rev. B 86, 245117 (2012);
P. V. Buividovich, Monte-Carlo study of quasiparticle dispersion relation in monolayer graphene,
PoS(Confinement X), 084 (2013).
[17] M. V. Ulybyshev, P. V. Buividovich, M. I. Katsnelson, and M. I. Polikarpov, Monte Carlo Study of the
Semimetal-Insulator Phase Transition in Monolayer Graphene with a Realistic Interelectron
Interaction Potential, Phys. Rev. Lett. 111, 056801 (2013).
[18] G. W. Semenoff, Condensed-matter Simulation of a Three-Dimensional Anomaly, Phys. Rev. Lett. 53,
2449 (1984).
[19] V. P. Gusynin, S. G. Sharapov, and J. P. Carbotte, AC conductivity of graphene: From tight-binding
model to 2+1-dimensional quantum electrodynamics, Int. J. Mod. Phys. B 21, 4611 (2007).
[20] C. Burden and A. N. Burkitt, Lattice Fermions in Odd Dimensions, Europhys. Lett. 3, 545 (1987).
[21] M. Göckeler, R. Horsley, P. Rakow, G. Schierholz, and R. Sommer, Scaling Laws, Renormalization
Group Flow and the Continuum Limit in Non-Compact Lattice QED, Nucl. Phys. B371, 713 (1992).
[22] M. Göckeler, R. Horsley, V. Linke, P. E. L. Rakow, G. Schierholz, and H. Stüben, Seeking the
equation of state of non-compact lattice QED, Nucl. Phys. B487, 313 (1997).
[23] I. F. Herbut, Quantum Critical Points with the Coulomb Interaction and the Dynamical Exponent:
When and Why z = 1, Phys. Rev. Lett. 87, 137004 (2001).
7

Download Report