Many-body instabilities and mass generation in slow Dirac materials

0.22 MB·English

by  Christopher Triola

Many-body instabilities and mass generation in slow Dirac materials
Christopher Triola,1,2 Jian-Xin Zhu,3,4 Albert Migliori,5 and Alexander V. Balatsky2,6
1Department of Physics, College of William and Mary, Williamsburg, Virginia 23187, USA
2Institute for Materials Science, Los Alamos National Laboratory, Los Alamos New Mexico 87545, USA
3Theoretical Division, Los Alamos National Laboratory, Los Alamos, New Mexico 87545, USA
4Center for Integrated Nanotechnology, Los Alamos National Laboratory, Los Alamos, New Mexico 87545, USA
5Glenn T. Seaborg Institute, Los Alamos National Laboratory, Los Alamos, New Mexico 87545, USA
5 6Nordic Institute for Theoretical Physics (NORDITA),
1
Roslagstullsbacken 23, S-106 91 Stockholm, Sweden
0
2 Some Kondo insulators are expected to possess topologically protected surface states with linear
g Diracspectrum,thetopologicalKondoinsulators. Becausethebulkstatesofthesesystemstypically
u have heavy effective electron masses, the surface states may exhibit extraordinarily small Fermi
A velocities that could force the effective fine structure constant of the surface states into the strong
couplingregime. Usingatight-bindingmodelwestudythemany-bodyinstabilitiesofthesesystems

3 andidentifyregions of parameterspacein whichthesystemexhibitsspindensitywave, andcharge
1 density waveorder.

]
l
l I. INTRODUCTION involving suspended graphene, for which α 2.2, seem
a ≈
toindicateagaplessstatetowithin0.1meVoftheDirac
h
- In many metals and semiconductors the behavior of point14. Therefore,thegroundstateofDiracmaterialsin
s the low energy electronic states can be understood in the strong coupling regime is not currently understood.
e
m terms of free quasiparticles with quadratic energy dis- For this reason we propose studying a class of materials
persion in momentum p2/2m∗, where m∗ is a system- with much smaller Fermi velocity than that of graphene
.
at dependent effective mass1. However,there are a number since this class of materials is likely to possess α ≫ αc
m of materials whose low energy electronic states are bet- andwouldbeabettercandidateforexperimentsprobing
ter described as massless Dirac fermions, including the the strong coupling regime in Dirac materials.
-
d superfluidphaseof3He,high-temperatured-wavesuper-
Thesurfaceofathree-dimensional(3D)topologicalin-
n conductors, graphene, and the surface states of topolog-
sulator (TI) hosts two-dimensional(2D) Dirac quasipar-
o ical insulators1–4. In these Dirac materials the kinetic
ticles similar to those found in graphene. Examples of
c
energy is proportional to the momentum vp, just like
[ experimentally verified 3D TIs include Bi Se , Bi Te ,
2 3 2 3
massless relativistic particles but with a speed v that
2 depends on the details of the system. For example, in ahnadlf oSfb2thTaet3,inalglroafphwehniech15,h16a.veHFoweremveirv,etlhoecriteieiss aronuogthhelyr
v graphene v 106m/s c/300.
8 ≈ ≈ class of topological insulators, the topological Kondo in-
The fact that quasiparticles obey the Dirac equation
5 sulators (TKI), in which the bulk states are formed by
insteadoftheSchr¨odingerequationcanaffectavarietyof
1 renormalizedf-electron levels which hybridize with con-
3 electronic properties, for example, the integer quantum ductionelectronstoformamillivolt-scalegapinthebulk
0 Halleffectandlocalization5. Anotherimportantwaythe spectrum17–20. The small gap in these materials com-
. Dirac nature of the quasiparticles manifests itself is in
1 bined with the large bulk effective mass imply that the
the effect of interactions. If the quasiparticles of a sys-
0 surface Fermi velocity could be quite small. Some ma-
5 tem obeythe Schr¨odingerequation,then the ratioofthe terials theoretically predicted to fall into this category
1 average interparticle Coulomb energy to the average ki- include SmB 21, YbB 22, and PuB 23. Furthermore,
: netic energy, r = E /E , is related to the density by 6 12 6
v s C K thereisagrowingbodyofexperimentalevidencedemon-
r n−1/d,6,7 where the constant of proportionality de-
Xi psen∝ds on characteristics of the material. In contrast to strating that SmB6 does in fact host metallic surface
states24–30.
r normalmetals, for Dirac materials this ratio is a charac-
a
teristicofthesystem,independentoftheelectrondensity, Previous work has explored the possibility of broken
given by α E /E = e2/(~ǫv). In this expression e symmetry states on the surface of a TKI employing a
C K
is the charg≡e of the electron, ǫ is the material’s dielec- continuum model31. In this paper we present a tight-
tric constant, ~ is the reduced Planck constant, and v is binding model to study the surface states of a TKI and
the speed of the Dirac particles. Much work has gone proceedtoinvestigatethepossibleorderedgroundstates
into the study ofthe phasediagramofgraphenewithre- for these systems within a mean field theory. From this
spect to this parameter α2,6–13 and the results indicate analysis,wefindregionsofparameterspaceforthemodel
that there is a critical value α such that if α < α the that admit spin density wave and charge density wave
c c
spectrumremainsgaplessandifα>α the systemflows solutions. For the case of strictly repulsive interactions
c
towardthestrongcouplingregimeandislikelytodevelop we find that these orderedsolutions lie within the region
agap6. Thusfar,perturbativeandnumericalresultssug- of parameter space correspondingto the strong coupling
gestthecriticalvalueisα 16,8,12,13,whileexperiments regime of Dirac materials (α>α 1).
c c
≈ ≈
2
II. THEORETICAL MODEL AND METHODS
To model the 2D surface states of a TKI, we consider
a Hamiltonian defined on a square lattice:
A
H = i ψ† zˆ Rˆ σ ψ
0 − 2 i,α,σ · ij × αβ j,β,σ
αX,β,σXhiji (cid:16) (cid:17) (1)
+ Γ ψ† σ ψ
ij i,α,σ αβ j,β,σ
αX,β,σXi,j FIG. 1. (Color online) Plots of the band structure given by
Eq. (3) along thediagonal of thesquareBrillouin zone using
where α and β are orbital indices, σ is a spin index, Rˆij four different values of the parameter A [1/4 (solid, black),
is the unit vector pointing from lattice site j to lattice 1/8 (dashed, red), 1/16 (dashed-dot, blue), and 1/32 (dot-
site i, and the matrix Γˆ is defined as ted, green)], in two different momentum ranges [(a) from
k = (−1/5a,−1/5a) to k = (1/5a,1/5a) and (b) from
4Γ ; i=j k = (−π/a,−π/a) to k = (π/a,π/a)]. All energies are in
Γ = Γ ; i,j nearest neighbors . (2) unitsof thebandwidth.
ij  −
0 ; otherwise

The term proportional to A leads to the formation of
four separate Dirac points in the Brillouin zone. The
Note that for A 4Γ the bandwidth is a constantsetby
term proportional to Γ acts as a momentum-dependent ≤
the model parameter Γ. In the analysis that follows we
mass term which gaps out all of the Dirac points except
restricttherangeofAtoA 4Γandpresentallenergies
the one at k = 0 allowing the model to represent the ≤
in units of the bandwidth w = 16Γ. We also present all
surface states of a strong TI32. The energy eigenvalues
distances in units of the lattice constant a.
associatedwiththisHamiltonianink spacearegivenby:
Toaccountforinteractionsweconsiderthe fullHamil-
ak ak tonian: H =H +H ,whereH takesontheexactform:
E± = 4Γ sin2 x +sin2 y 0 I I
k
± 2 2
(cid:20) (cid:21)
×vuu1+(cid:18)4AΓ(cid:19)2 ssiinn22aak2kxx++ssiinn22aak2kyy 2. (3) HI = V20 Xσ,σ′Xi6=j |er−i|−ri−rjrj||2/λ+d2ψi†,f,σψi,f,σψj†,f,σ′ψj,f,σ′
ut h i −U ψi†,f,↑ψpi,f,↑ψi†,f,↓ψi,f,↓
Expandingthisdispersionforsmall k alongthekx direc- Xi
tion we find: (5)
whereV controlsthestrengthofthelong-rangeCoulomb
0
(3Γ2 A2)a3 interaction between f electrons and U is introduced as
E± aAk+ − k3 . (4)
k ≈± 6A an on-site interaction between f electrons.
(cid:18) (cid:19)
In our calculations, we replace the exact interaction
Thus, we can see that to first order in k the dispersion
term H with the mean field Hamiltonian:
matches the Dirac dispersion with Fermi velocity given I
by aA/~. In Fig. 1 we plot the full dispersion from Eq.
(3)fordifferentrangesofk todemonstratetheDiracdis- e−|ri−rj|/λ
persionforafewdifferentvaluesofthe Fermivelocity. It HMF =V n ψ† ψ
shows that near the Dirac point the parameter A con- I 0Xσ,σ′Xi6=j |ri−rj |2 +d2h i,f,σi j,f,σ′ j,f,σ′
trols the Fermi velocity; however, for A Γ we can see
≪ U n p ψ† ψ + n ψ† ψ
that the cubic term in Eq. (4) begins to dominate and − h i,f,↑i i,f,↓ i,f,↓ h i,f,↓i i,f,↑ i,f,↑
the dispersion away from the Dirac point becomes no- Xi (cid:16) (cid:17)
ticeably less linear. Since we are most interested in the + ∆ ψ† ψ† +∆∗ψ ψ +E .
i i,f,↑ j,f,↓ i i,f,↓ j,f,↑ 0
regime in which the model best describes a Dirac ma-
Xi (cid:16) (cid:17)
terial, in this paper we focus on the case in which the (6)
chemical potential is close to the Dirac point. Tocapturetheon-siteCoulombrepulsionbetweenf elec-
In terms of the model parameters the bandwidth is trons we can set U = V0/d. If we wish to include an
−
given by: attractive on-site interaction we set U > 0. In the ab-
sence of a large electron-phonon coupling this attractive
16Γ ; A 4Γ interaction could be engineered by the adsorption of a
≤
w = 4A2 finite density ofnonmagneticmolecules asdiscussedpre-
 ; A>4Γ
 √2A2+16Γ2 viously in the context of topological insulators33. We

3
may choose to write this in a more compact notation as: where f(ǫ) = 1 is the Fermi-Dirac distribution
eǫ/kBT+1
functionattemperatureT andk istheBoltzmanncon-
B
HMF = W ψ† ψ stant. Given an initial set of model parameters and
I iασ,jβσ′ j,β,σ′ j,β,σ′
α,β,σ,σ′ i,j a temperature, Eqs. (7) and (8) allow us to solve for
X X
the density profile and superconducting order parameter
+ ∆ ψ† ψ† +∆∗ψ ψ +E
i i,f,↑ j,f,↓ i i,f,↓ j,f,↑ 0 ∆ self-consistently. In the next section we discuss our
Xi (cid:16) (cid:17) progress toward solving these equations.
In some cases multiple solutions for the same model
where
parameters may be found. In this case it is useful to
V e−|ri−rj|/λ n δ δ ; i=j compare the free energy associated with each of the so-
Wiασ,jβσ′ = 0√|ri−rj|2+d2h i,f,σi αβ αf 6 lutions, given by F =kBT lnZ, where Z is the partition
( U ni,f,σ δαβδαf(1 δσσ′) ; i=j function. The true ground state of the system will be
− h i −
given by the solution with the lowest free energy.
and
∆ U ψ ψ .
i ≡ h i,f,↑ j,f,↓i III. NUMERICAL RESULTS AND DISCUSSION
Equippedwiththismean-fieldHamiltonianweperform
While it is straightforward to numerically solve Eqs.
a Bogoliubov transformation:
(7) and (8) for a finite system, we can make the compu-
ψ = γ u γ† v∗ tation more efficient by using the supercell technique34
i,α,↑ n↑ i,α,n,↑− n↓ i,α,n,↑ asdescribedintheappendix. Forasystemwitha10 10
n ×
X real space unit cell and an 8 8 supercell we solved for
ψi,α,↓ = γn↓ui,α,n,↓+γn†↑vi∗,α,n,↓, self-consistent solutions to E×qs. (7) and (8). For the
n followingresultswe focusedonthe caseofnearestneigh-
X
borCoulombinteractionsonly andthe zerotemperature
where γ† (γ ) creates (annihilates) an eigenstate of
nσ nσ limit. InEq(6)weusedascreeninglengthofλ=1anda
the mean-fieldHamiltonianH. It can be shownthat the
lattice cutoff ofd=1. We consideredtwo limiting cases:
coefficients u and v satisfy the following equations:
thecaseofarepulsiveon-siteinteraction(U = V ),and
0
−
the case of an attractive on-site interaction that scales
ǫn,↑ui,α,n,↑ = Hiα↑,jβ↑uj,β,n,↑+∆ivi,f,n,↓ with the Coulomb interaction (U =V0).
j,β Startingfrominitialseedsthatpossessedantiferromag-
X
ǫ v = H∗ v +∆∗u netic,ferromagnetic,checkerboardandstripechargeden-
n,↑ i,α,n,↓ − iα↓,jβ↓ j,β,n,↓ i i,f,n,↑ sity wave (CDW) order in addition to random seeds we
j,β
X (7) found self-consistent solutions for Eqs. (7) and (8) us-
ǫn,↓ui,α,n,↓ = Hiα↓,jβ↓uj,β,n,↓+∆ivi,f,n,↑ ing a convergence criterion of 10−3. Some of the self-
j,β consistentsolutionsthatemergedfromthedifferentseeds
X
ǫ v = H∗ v +∆∗u for the same model parameters differed from each other.
n,↓ i,α,n,↑ − iα↑,jβ↑ j,β,n,↑ i i,f,n,↓ In these cases the one with the lowest free energy was
j,β
X taken to be the solution. In Fig. 2 we show the regions
of parameter space for which we found solutions in the
where H H(0) + W and ǫ are
iασ,jβσ′ ≡ iασ,jβσ′ iασ,jβσ′ n,σ case of repulsive on-site interactions while in Fig. 4 we
eigenvalues of H.
show the regions of parameter space for which we found
Giventhesolutionstotheseequationswecanwritethe
solutions in the case with attractive on-site interactions.
mean fields as
First, we consider the case of on-site repulsion (U =
V ), Fig. 2. Note that the general trend is consistent
hni,α,↑i= |ui,α,n,↑ |2 f(ǫn,↑) −wit0h our expectations for Dirac materials. In the region
n
X of strong coupling, α = V /A > α , we find Coulomb-
0 c
+ vi,α,n,↓ 2 (1 f(ǫn,↓)) drivenorderedstates, while in the weak coupling region,
| | −
Xn V0/A < αc, a paramagnetic (PM) normal metallic state
n = u 2 f(ǫ ) exists. These results are consistent with the established
i,α,↓ i,α,n,↓ n,↓
h i | | value of α 1. However, it appears that there is a
Xn (8) critical valcue≈of the coupling, V w/3, for this model
+ v 2 (1 f(ǫ )) c ≈
| i,α,n,↑ | − n,↑ below which the solution is trivial. This is in contrast
Xn to the case of a Dirac continuum model in which the
∆i =U vi∗,f,n,↓ui,f,n,↑(1−f(ǫn,↑)) only parameter governing the Coulomb interaction is α.
n Thisdifferencecanbeattributedtothefactthatforvery
X
U v∗ u f(ǫ ) small values of A the band structure appears less linear
− i,f,n,↑ i,f,n,↓ n,↓ and eventually the cubic term becomes more important,
n
X
4
with (π,π) wave vectors as shown in the sample plot in
Fig. 3. Intermediate states were also observed but these
appear to be higher energy excitations. In the SDW re-
gion the boundary for the phase along the V axis, at
0
approximately one third of the bandwidth, defines the
criticalcoupling, V . We find that above another critical
c
value of V CDW orderbegins to coexist with the SDW.
0
In the coexistence region for some model parameters we
were able to find solutions with exclusively CDW order
but we lacked the resolution to see if these solutions in-
dicatedthe existence ofanadditionalregionofthe plane
in which CDW is truly favored over SDW order, further
calculations will be needed to answer this question.
Next we turn our attention to the case in which we
include an on-site attraction (U = V ), as shown in Fig.
0
4. In this case we find two regions: a region with CDW
andaPMregion. Again,thesedensitymodulationswere
associated with a wave vector of (π,π). The region that
favors CDW begins at V w/3 and covers the rest of
0
≈
the plane. It is interesting to note that the CDW or-
der appears for V > w/3 which is the same as V for
0 c
FIG.2. (Coloronline)Plotofthephasesfortheself-consistent the case with repulsive on-site interactions. It should be
solutions found in different regions of the A,V0 plane. Note notedthatsome ofthe self-consistentsolutionswe found
that the region below the line A = V0 appears to favor the near the transition region V w/3 seemed to possess
formation of nontrivial order, consistent with αc ≈ 1. The 0 ≈
a small superconducting order parameter; however, this
regionenclosedbythereddashedlinefavorstheformationof
order parameter was usually just below the convergence
spin density wave order while in the region enclosed by the
criterion (even when the convergence criterion was low-
black solid line we find both spin density wave and charge
density wave solutions. Outside of these regions the solution ered to 10−7). We attribute the absence of a supercon-
is paramagnetic (PM). ducting region to the fact that we restricted ourselves
to the case of half filling in which there was no density
of states to allow for superconducting pairing. A more
detailed study of the region near V w/3 may be in-
0
aswecanseefromEq. (4). Itisreasonabletoexpectthat teresting for future work studying this≈model away from
realmaterials which host slow Dirac states will typically half filling.
havesimilarbehaviorsince the bands forthese materials
NoteinFig. 4theabsenceofanyregionswithmagnetic
areexpectedtodevelopnonzerocurvatureawayfromthe
order,in contrastto Fig. 2 in which both AFM and FM
Dirac point18,19,24,29.
orderwerefound. Thiscanbe accountedforbyaheuris-
tic argument based on Eq. (6). Notice that the spin-
dependenttermsinthemeanfieldHamiltonianaregiven
by U n ψ† ψ + n ψ† ψ ,thus
− i h i,f,↑i i,f,↓ i,f,↓ h i,f,↓i i,f,↑ i,f,↑
the expect(cid:16)ation value of the contribution to the(cid:17) total
P
energy will be 2U n n . For U > 0 we
− ih i,f,↑ih i,f,↓i
can see that the energy can be minimized if the sum
P
n n takes on its maximum possible value.
ih i,f,↑ih i,f,↓i
Each term of this sum has a maximum value when
P
n = n =1/2. Therefore the minimum energy
i,f,↑ i,f,↓
h i h i
can be expected to be achieved in a state with no mag-
netic order. However, for U < 0 the minimum energy is
FIG.3. (Coloronline)Plotofthedensitymodulationsovera achievedforaminimumvalueof n n ,which
10×10realspaceunitcell,asobservedintheSDWandCDW can allow the system to minimize iiths ie,nf,e↑righyit,fh,↓roiugh an
regions shown in Figs. 2 and 4. on-site spin polarization. P
TakingacloserlookatFig. 2wecanseethatthereare
IV. CONCLUSIONS
threedistinctregionsoftheV ,Aplane: aregionfavoring
0
spin density wave (SDW) order, a region in which SDW
andCDWcoexist,andaregioninwhichthesolutionwas In conclusion, we presented a model for studying the
PM. Both SDW and CDW modulations were associated surface states of a class of topological Kondo insulators
5
Science at Los Alamos National Laboratory. This re-
search used resources of the National Energy Research
Scientific Computing Center (NERSC), a DOE Office of
Science User Facility supported by the Office of Science
of the U.S. Department of Energy under Contract No.
DE-AC02-05CH11231.
Appendix A: Supercell Technique
ThesystemofequationsgivenbyEqs. (7)and(8)can
be solved for finite systems by simple matrix diagonal-
ization. However, the matrix that must be diagonalized
is 8N 8N, where N is the number of lattice sites and
×
8 = 2(spins) 2(orbitals) 2(electron-hole). We can
× ×
see that for 6400 sites this would involve diagonalizing
a 51200 51200 matrix which is not terribly practical.
×
Using the supercell technique we can decrease the size
of the matrix that needs to be diagonalizedsignificantly.
Inthe frameworkofthe supercelltechniquewe recognize
that, due to the periodicity of the system, the solutions
FfoIuGn.d4i.ndPilffoetreonfttrheegiponhsasoefstfhoerAth,Ve0spellaf-nceonfosirstaetnttrascotliuvetioonns- tuhriis,αw,ne,σwarnitdevri,α,n,σ, are Bloch waves. To account for
site interaction. In the region enclosed to the right of the
solid black line the self-consistent solutions possessed charge uri,α,n,σ =ei∗ri·kuk,ri,α,n,σ (A1)
dpeanrasmityagwnaevtieco(PrdMer);. outside of this region the solution was vri,α,n,σ =ei∗ri·kvk,ri,α,n,σ,
where k is the crystal momentum. After this transfor-
mation Eqs. (7) and (8) become:
and explored the dependence of the band structure on ǫk,n,↑uk,i,α,n,↑ = Hiα↑,jβ↑;kuk,j,β,n,↑+∆ivk,i,f,n,↓
the model parameters, identifying the parameters which j,β
X
determine the Fermi velocity at the Dirac point. We
then added interactions to this model, accounting for ǫk,n,↑vk,i,α,n,↓ =− Hi∗α↓,jβ↓;kvk,j,β,n,↓+∆∗iuk,i,f,n,↑
j,β
both Coulomb interactions as well as the possibility of X
an on-site attractive interaction. Using mean-field the- ǫk,n,↓uk,i,α,n,↓ = Hiα↓,jβ↓;kuk,j,β,n,↓+∆ivk,i,f,n,↑
ory, at zero temperature, we found self-consistent solu- j,β
X
tionsfordifferentmodelparameters,investigatingthere-
lationshipbetweenthe Fermivelocityatthe Diracpoint, ǫk,n,↓vk,i,α,n,↑ =− Hi∗α↑,jβ↑;kvk,j,β,n,↑+∆∗iuk,i,f,n,↓
j,β
the strength of the interactions, and the nature of the X
(A2)
self-consistent solutions. For the case with on-site repul-
and
sion we identified three regions of parameter space with
1
different Fermi velocity and coupling strength: a region ni,α,↑ = uk,i,α,n,↑ 2 f(ǫk,n,↑)
which exclusively favored spin density wave order, a re- h i Mxy n,k | |
gion of coexisting spin density wave and charge density X
1
wave order, and a paramagnetic normal metallic region. + vk,i,α,n,↓ 2 (1 f(ǫk,n,↓))
M | | −
We also identified a critical value of the Coulomb in- xy n,k
X
teraction strength Vc w/3 below which the solutions 1
were normal metallic.≈When we considered the case of ni,α,↓ = uk,i,α,n,↓ 2 f(ǫk,n,↓)
h i M | |
an attractive on-site interaction we found that the so- xy n,k
X
lutions possessed charge density wave order above this 1
same critical Coulomb interaction strength. + M |vk,i,α,n,↑ |2 (1−f(ǫk,n,↑))
xy n,k
Acknowledgements: We are grateful to Yuan-Yen Tai, X
U
TLuowo,fiqEAnrhimcoedR,oDsasiv,idaAndberZgaeclh,aYroynaFtiasnkDfuorbi,uYseofnuglkadnisg- ∆i = Mxy n,kvk∗,i,f,n,↓uk,i,f,n,↑(1−f(ǫk,n,↑))
cussions. This work was supported by US DOE BES X
U
E304. The work of A.V.B. was also supported by ERC − M vk∗,i,f,n,↑uk,i,f,n,↓f(ǫk,n,↓),
DM321031andKAW.C.T.acknowledgesthe hospitality xy n,k
X
of the Seaborg Institute and the Institute for Materials (A3)
6
where k = 2π nx, ny where n = 1,2,...,M and Now, a system composed of 6400 sites can be studied
Mxya Nx Ny x x by diagonalizing a 10 10 real space system using an
ny =1,2,...,My,M(cid:16) x andM(cid:17)y arethenumberofunitcells ×
8 8 supercell. This means we only need to diagonal-
in the x and y direction, respectively, Mxy =MxMy, Nx ×
ize a 800 800 matrix instead of 51200 51200. More-
andNy arethenumberoflatticesitesperunitcellinthe over, this×diagonalization is performed f×or each k inde-
x and y direction, respectively, and we define
pendently and thus the procedure may be easily paral-
Hiασ,jβσ′;k ≡ eik·(rj+Rj−ri)Hriασ,(rj+Rj)βσ′. lelized to further improve performance.
R
Xj
1 T. Wehling, A. Black-Schaffer, and A. Balatsky, Adv. 19 M.Dzero,K.Sun,V.Galitski,andP.Coleman,Phys.Rev.
Phys. 63, 1 (2014). Lett. 104, 106408 (2010).
2 H. P. Dahal, T. O. Wehling, K. S. Bedell, J.-X. Zhu, and 20 B. Roy, J. D. Sau, M. Dzero, and V. Galitski, Phys. Rev.
A. Balatsky, Physica B 405, 2241 (2010). B 90, 155314 (2014).
3 L. Fu, C. L. Kane, and E. J. Mele, Phys. Rev. Lett. 98, 21 V. Alexandrov, M. Dzero, and P. Coleman, Phys. Rev.
106803 (2007). Lett. 111, 226403 (2013).
4 L. Fu and C. L. Kane, Phys.Rev.B 76, 045302 (2007). 22 H. Weng, J. Zhao, Z. Wang, Z. Fang, and X. Dai, Phys.
5 A. C. Neto, F. Guinea, N. Peres, K. S. Novoselov, and Rev.Lett. 112, 016403 (2014).
A. K.Geim, Rev.Mod. Phys.81, 109 (2009). 23 X.Deng, K.Haule, and G. Kotliar, Phys.Rev.Lett.111,
6 V. N. Kotov, B. Uchoa, V. M. Pereira, F. Guinea, and 176404 (2013).
A. C. Neto, Rev.Mod. Phys.84, 1067 (2012). 24 M. Neupane, N. Alidoust, S.-Y. Xu, T. Kondo, Y. Ishida,
7 H. P. Dahal, Y. N. Joglekar, K. S. Bedell, and A. V. Bal- D. Kim, C. Liu, I. Belopolski, Y. Jo, T.-R. Chang, et al.,
atsky,Phys. Rev.B 74, 233405 (2006). Nat. Commun. 4 (2013).
8 J.E.DrutandT.A.L¨ahde,Phys.Rev.Lett.102,026802 25 J. Jiang, S. Li, T. Zhang, Z. Sun, F. Chen, Z. Ye, M. Xu,
(2009). Q.Ge,S.Tan,X.Niu,etal.,Nat.Commun.4,3010(2013).
9 D. Khveshchenko, J. Phys.: Condens. Matter 21, 075303 26 D. Kim, S. Thomas, T. Grant, J. Botimer, Z. Fisk, and
(2009). J. Xia, Sci. Rep.3, 3150 (2013).
10 S.Ryu,C.Mudry,C.-Y.Hou,andC.Chamon,Phys.Rev. 27 N. Xu, X. Shi, P. K. Biswas, C. E. Matt, R. S. Dhaka,
B 80, 205319 (2009). Y. Huang, N. C. Plumb, M. Radovi, J. H. Dil, E. Pom-
11 E. Gorbar, V. Gusynin, V. Miransky, and I. Shovkovy, jakushina, et al., Phys. Rev.B 88, 121102 (2013).
Phys. Rev.B 66, 045108 (2002). 28 S. Wolgast, C¸. Kurdak, K. Sun, J. Allen, D.-J. Kim, and
12 O. Gamayun, E. Gorbar, and V. Gusynin, Phys. Rev. B Z. Fisk, Phys. Rev.B 88, 180405 (2013).
81, 075429 (2010). 29 N. Xu, P. K. Biswas, J. H. Dil, R. S. Dhaka, G. Landolt,
13 O. Gamayun, E. Gorbar, and V. Gusynin, Phys. Rev. B S. Muff, C. E. Matt, X. Shi, N. C. Plumb, M. Radovi,
80, 165429 (2009). et al., Nat.Commun. 5, 4566 (2014).
14 D. Elias, R. Gorbachev, A. Mayorov, S. Morozov, 30 D.-J.Kim,J.Xia,andZ.Fisk,Nat.Mater.13,466(2014).
A. Zhukov, P. Blake, L. Ponomarenko, I. Grigorieva, 31 B.Roy,J.Hofmann,V.Stanev,J.D.Sau,andV.Galitski,
K. Novoselov, F. Guinea, et al., Nat. Phys.7, 701 (2011). arXiv:1410.1868 (2014).
15 H. Zhang, C.-X. Liu, X.-L.Qi, X. Dai, Z. Fang, and S.-C. 32 O.Vafek andA.Vishwanath,Ann.Rev.Condens.Matter
Zhang, Nature Physics5, 438 (2009). Phys. 5, 83 (2014).
16 J. Zhang, C. Triola, and E. Rossi, Phys. Rev. Lett. 112, 33 J.-H. She and A. V. Balatsky, Phys. Rev. B 90, 104517
096802 (2014). (2014).
17 P. Coleman, Handbook of magnetism and advanced mag- 34 J.-X. Zhu, B. Friedman, and C. Ting, Phys. Rev. B 59,
netic materials, vol 1, 95 (2007). 3353 (1999).
18 M.Dzero,K.Sun,P.Coleman,andV.Galitski,Phys.Rev.
B 85, 045130 (2012).