Electron states in the field of charged impurities in two-dimensional Dirac systems
We review the theoretical and experimental results connected with the electron states in two-dimensional Dirac systems paying a special attention to the atomic collapse in graphene. Two-electron bound states of a Coulomb impurity are considered too. A rather subtle role of a magnetic field in the su...
Gespeichert in:
| Datum: | 2018 |
|---|---|
| Hauptverfasser: | , , |
| Format: | Artikel |
| Sprache: | English |
| Veröffentlicht: |
Фізико-технічний інститут низьких температур ім. Б.І. Вєркіна НАН України
2018
|
| Schriftenreihe: | Физика низких температур |
| Schlagworte: | |
| Online Zugang: | https://nasplib.isofts.kiev.ua/handle/123456789/176115 |
| Tags: |
Tag hinzufügen
Keine Tags, Fügen Sie den ersten Tag hinzu!
|
| Назва журналу: | Digital Library of Periodicals of National Academy of Sciences of Ukraine |
| Zitieren: | Electron states in the field of charged impurities in two-dimensional Dirac systems / E.V. Gorbar, V.P. Gusynin, O.O. Sobol // Физика низких температур. — 2018. — Т. 44, № 5. — С. 491-524. — Бібліогр.: 142 назв. — англ. |
Institution
Digital Library of Periodicals of National Academy of Sciences of Ukraine| id |
nasplib_isofts_kiev_ua-123456789-176115 |
|---|---|
| record_format |
dspace |
| spelling |
nasplib_isofts_kiev_ua-123456789-1761152025-02-10T01:52:07Z Electron states in the field of charged impurities in two-dimensional Dirac systems Gorbar, E.V. Gusynin, V.P. Sobol, O.O. Обзор We review the theoretical and experimental results connected with the electron states in two-dimensional Dirac systems paying a special attention to the atomic collapse in graphene. Two-electron bound states of a Coulomb impurity are considered too. A rather subtle role of a magnetic field in the supercritical charge problem in graphene is discussed. The electron states in the field of two equally charged impurities are studied and the conditions for supercritical instability to occur are determined. It is shown that the supercriticality of novel type is realized in gapped graphene with two unlikely charged impurities. For sufficiently large charges of impurities, it is found that the wave function of the occupied electron bound state of the highest energy changes its localization from the negatively charged impurity to the positively charged one as the distance between the impurities increases. The specifics of the atomic collapse in bilayer graphene is considered and it is shown that the atomic collapse in this material is not related to the phenomenon of the fall-to-center. The work of E.V.G. and V.P.G. is partially supported by the National Academy of Sciences of Ukraine (project 0116U003191) and by its Program of Fundamental Research of the Department of Physics and Astronomy (project No. 0117U000240). V.P.G. acknowledges the support of the RISE Project CoExAN GA644076. 2018 Article Electron states in the field of charged impurities in two-dimensional Dirac systems / E.V. Gorbar, V.P. Gusynin, O.O. Sobol // Физика низких температур. — 2018. — Т. 44, № 5. — С. 491-524. — Бібліогр.: 142 назв. — англ. 0132-6414 PACS: 71.70.Di, 73.22.Pr, 81.05.ue https://nasplib.isofts.kiev.ua/handle/123456789/176115 en Физика низких температур application/pdf Фізико-технічний інститут низьких температур ім. Б.І. Вєркіна НАН України |
| institution |
Digital Library of Periodicals of National Academy of Sciences of Ukraine |
| collection |
DSpace DC |
| language |
English |
| topic |
Обзор Обзор |
| spellingShingle |
Обзор Обзор Gorbar, E.V. Gusynin, V.P. Sobol, O.O. Electron states in the field of charged impurities in two-dimensional Dirac systems Физика низких температур |
| description |
We review the theoretical and experimental results connected with the electron states in two-dimensional Dirac systems paying a special attention to the atomic collapse in graphene. Two-electron bound states of a Coulomb impurity are considered too. A rather subtle role of a magnetic field in the supercritical charge problem in graphene is discussed. The electron states in the field of two equally charged impurities are studied and the conditions for supercritical instability to occur are determined. It is shown that the supercriticality of novel type is realized in gapped graphene with two unlikely charged impurities. For sufficiently large charges of impurities, it is found that the wave function of the occupied electron bound state of the highest energy changes its localization from the negatively charged impurity to the positively charged one as the distance between the impurities increases. The specifics of the atomic collapse in bilayer graphene is considered and it is shown that the atomic collapse in this material is not related to the phenomenon of the fall-to-center. |
| format |
Article |
| author |
Gorbar, E.V. Gusynin, V.P. Sobol, O.O. |
| author_facet |
Gorbar, E.V. Gusynin, V.P. Sobol, O.O. |
| author_sort |
Gorbar, E.V. |
| title |
Electron states in the field of charged impurities in two-dimensional Dirac systems |
| title_short |
Electron states in the field of charged impurities in two-dimensional Dirac systems |
| title_full |
Electron states in the field of charged impurities in two-dimensional Dirac systems |
| title_fullStr |
Electron states in the field of charged impurities in two-dimensional Dirac systems |
| title_full_unstemmed |
Electron states in the field of charged impurities in two-dimensional Dirac systems |
| title_sort |
electron states in the field of charged impurities in two-dimensional dirac systems |
| publisher |
Фізико-технічний інститут низьких температур ім. Б.І. Вєркіна НАН України |
| publishDate |
2018 |
| topic_facet |
Обзор |
| url |
https://nasplib.isofts.kiev.ua/handle/123456789/176115 |
| citation_txt |
Electron states in the field of charged impurities in two-dimensional Dirac systems / E.V. Gorbar, V.P. Gusynin, O.O. Sobol // Физика низких температур. — 2018. — Т. 44, № 5. — С. 491-524. — Бібліогр.: 142 назв. — англ. |
| series |
Физика низких температур |
| work_keys_str_mv |
AT gorbarev electronstatesinthefieldofchargedimpuritiesintwodimensionaldiracsystems AT gusyninvp electronstatesinthefieldofchargedimpuritiesintwodimensionaldiracsystems AT soboloo electronstatesinthefieldofchargedimpuritiesintwodimensionaldiracsystems |
| first_indexed |
2025-12-02T14:24:09Z |
| last_indexed |
2025-12-02T14:24:09Z |
| _version_ |
1850406809714556928 |
| fulltext |
Low Temperature Physics/Fizika Nizkikh Temperatur, 2018, v. 44, No. 5, pp. 491–524
Electron states in the field of charged impurities
in two-dimensional Dirac systems
(Review Article)
E.V. Gorbar1,2, V.P. Gusynin2, and O.O. Sobol1
1Department of Physics, Taras Shevchenko National University of Kiev, Kiev 03680, Ukraine
E-mail: vgusynin@bitp.kiev.ua
2Bogolyubov Institute for Theoretical Physics, Kiev, 03680, Ukraine
Received February 18, 2018, published online March 27, 2018
We review the theoretical and experimental results connected with the electron states in two-dimensional Di-
rac systems paying a special attention to the atomic collapse in graphene. Two-electron bound states
of a Coulomb impurity are considered too. A rather subtle role of a magnetic field in the supercritical charge
problem in graphene is discussed. The electron states in the field of two equally charged impurities are studied
and the conditions for supercritical instability to occur are determined. It is shown that the supercriticality
of novel type is realized in gapped graphene with two unlikely charged impurities. For sufficiently large charges
of impurities, it is found that the wave function of the occupied electron bound state of the highest energy changes
its localization from the negatively charged impurity to the positively charged one as the distance between
the impurities increases. The specifics of the atomic collapse in bilayer graphene is considered and it is shown
that the atomic collapse in this material is not related to the phenomenon of the fall-to-center.
PACS: 71.70.Di Landau levels;
73.22.Pr Electronic structure of graphene;
81.05.ue Graphene.
Keywords: graphene, Dirac mass gap, atomic collapse, two-center problem.
Contents
1. Introduction .......................................................................................................................................... 492
2. Atomic collapse in monolayer graphene .............................................................................................. 493
2.1. Resonance states in gapless graphene in quasiclassical approach ................................................ 493
2.2. Supercritical instability in graphene ............................................................................................. 494
2.2.1. Gapped graphene, subcritical regime ................................................................................. 494
2.2.2. Gapped graphene, supercritical regime .............................................................................. 496
2.2.3. Gapless graphene ............................................................................................................... 497
2.3. Experimental observation of the atomic collapse in graphene ..................................................... 498
3. Supercritical instability in a magnetic field .......................................................................................... 499
3.1. The Coulomb center in a magnetic field....................................................................................... 500
3.2. Tuning the screening of charged impurity with chemical potential .............................................. 503
3.3. Screening charged impurities and lifting the orbital degeneracy in graphene by populating
Landau levels ............................................................................................................................... 505
4. Two-electron bound states near a Coulomb impurity ........................................................................... 506
5. Two-center problem ............................................................................................................................. 508
5.1. LCAO approach for symmetric two-center problem .................................................................... 508
5.2. Variational method ....................................................................................................................... 509
5.3. Quasistationary states ................................................................................................................... 512
6. The dipole problem and migration of the wave function ...................................................................... 512
6.1. Point electric dipole in graphene .................................................................................................. 513
6.2. Migration of the wave function in the finite dipole potential ....................................................... 514
7. Coulomb center problem in bilayer graphene ...................................................................................... 517
7.1. Two-band model .......................................................................................................................... 518
7.2. Numerical results ......................................................................................................................... 519
8. Summary .............................................................................................................................................. 521
References ................................................................................................................................................ 522
© E.V. Gorbar, V.P. Gusynin, and O.O. Sobol, 2018
E.V. Gorbar, V.P. Gusynin, and O.O. Sobol
1. Introduction
The phenomenon of the fall-to-center is deeply rooted
in the history of physics. The Rutherford’s discovery of the
planetary model of the atom immediately brought to the
light the problem of the stability of the atom. Indeed, clas-
sically, the electron rotating around the nucleus should emit
electromagnetic radiation, lose its energy, and fall to the
nucleus. We know that the atoms are stable and the atomic
collapse is avoided due to the uncertainty principle of
quantum mechanics. While the Coulomb interaction scales
like /Ze r− , where r is the distance to the nucleus and Ze
is its charge, the positive electron kinetic energy diverges
more strongly 2 2 2/ (2 ) / (2 )e em m rp as 0r → . There-
fore, the fall to the nucleus is energetically forbidden.
This qualitative argument shows that the fall-to-center
may still be possible in quantum mechanics for more sin-
gular potentials ( ) 1/ nV r r with 2n ≥ . In fact, the Schrö-
dinger equation with the potential 2( ) = /V r r−β provides
the canonical textbook example of the fall-to-center in quan-
tum mechanics [1], which takes place for 2> / (8 )emβ
when the energy spectrum is not bounded from below. If
the interaction potential is regularized at some distance 0r ,
then the electron wave function of the ground state is local-
ized in the region of the radius 0r which shrinks to the
origin as 0 0r → .
Still physically as | |p attains the value of order em c,
where c is the speed of light, the relativistic effects become
relevant. Since the kinetic term in the Dirac equation de-
pends linearly on momentum, the kinetic energy of the elec-
tron in the relativistic regime scales like / r as 0r → .
This means that already the Coulomb interaction could lead
to the atomic collapse. In quantum electrodynamics (QED)
for the regularized Coulomb potential, the atomic collapse
takes place for 170Z [2–4] when the lowest energy elec-
tron bound state dives into the lower continuum transform-
ing into a narrow resonance. This leads to the spontaneous
creation of electron-positron pairs with the electrons screen-
ing the positively charged nucleus and the positrons emit-
ted to infinity. Since supercritically charged nuclei are not
encountered in nature, this phenomenon was never observ-
ed in QED. It was suggested in the 70-ties [3,5–7] that the
supercritical instability in QED can nevertheless be exper-
imentally tested in a collision of two heavy nuclei. Alt-
hough subsequent experiments confirmed the existence of
supercritical fields in collisions of very heavy nuclei and
the gross features of positron emission [4], the analysis of
the supercritical regime turned out to be a difficult problem
mainly due to the transient nature of supercritical fields
generated during collisions.
It is an interesting question whether the supercritical in-
stability could be observed in the condensed matter sys-
tems. The first natural place to look is the narrow gap sem-
iconductors whose conductance and valence bands are
separated by a small gap. There exist also condensed mat-
ter systems with the relativistic-like energy spectrum of
quasiparticles. Bismuth, whose quasiparticles are described
by the massive Dirac equation, provides the historically first
example of such a system (for a review, see Refs. 8, 9).
Long time ago Herring argued [10] that the conductance
and valence bands in solids could, in general, meet at dis-
crete touching points. Remarkably, the energy dispersion
in the vicinity of these bands touching points is linear and
resembles the Weyl equation. The recently discovered Dirac
and Weyl semimetals whose itinerant electrons are descri-
bed by the 3D Dirac and Weyl equations, respectively, ex-
perimentally realize the Herring’s prediction (for a review,
see, e.g., Ref. 11). However, the corresponding materials are
characterized by the large dielectric constants. The small
value of the effective coupling constant makes it practically
impossible to realize the supercritical instability in these
materials.
The situation is different in graphene whose effective
coupling constant 2= e /( ) 2.2g Fα ≈v , where /300F c≈v
is the Fermi velocity, exceeds unity. This drastically de-
creases the value of the critical charge in graphene [12–15].
Although, according to the theory, the supercritical insta-
bility should be easily realized for charged impurities in
graphene, its experimental observation remained elusive
until recently. The problem is that it is difficult to produce
highly charged impurities because of their fast recombina-
tion. Still one can reach the supercritical regime by collect-
ing a large enough number of charged impurities in a cer-
tain region of graphene. Such an approach was recently
successfully realized [16] by using the tip of a scanning
tunneling microscope in order to create clusters of charged
calcium dimers.
In addition, the external charge in the realistic experi-
mental set-up should be smeared over a finite region of the
graphene plane because, otherwise, the Dirac equation is
no longer applicable and other nearest σ-bands should be
included in the analysis [15]. Thus, the potential of charged
impurities should be necessarily regularized at small dis-
tances in order that the continuum problem be well posed
physically. For instance, the charged impurities displaced
from the graphene plane provide such a natural regulariza-
tion and help to avoid the reconstruction of the band spec-
trum which takes place if they are placed directly into the
graphene plane or a disorder is present [17–20].
An interesting aspect of the electron physics in graphene
is its two-dimensional character. Therefore, the supercriti-
cal instability in the field of a charged impurity in graphene
is, in fact, the atomic collapse in a Flatland. Of course, this
does not mean that the theory governing the electron-elec-
tron interactions in graphene is QED in (2+1) dimensions.
Although the electrons are confined in the plane of graphene,
the electromagnetic force lines spread beyond the graphene’s
plane resulting in the standard Coulomb interaction poten-
tial 2( ) = e /CV r r . The crucial advantage of graphene com-
pared to QED is its experimental accessibility where atom-
492 Low Temperature Physics/Fizika Nizkikh Temperatur, 2018, v. 44, No. 5
Electron states in the field of charged impurities in two-dimensional Dirac systems
ic collapse can be investigated in table-top experiments
varying such parameters as doping and gate voltage.
The supercritical charge instability is closely related to
the excitonic instability in graphene in the strong coupling
regime > 1cα α (see, Refs. 21–23) and possible gap
opening, which may transform graphene into an insulator
[24–33]. Indeed, the excitonic instability can be viewed as
a many-body analog of the supercritical instability in the
field of a charged impurity and the critical coupling cα is
an analog of the critical coupling constant cZ α in the prob-
lem of the Coulomb center. In the strong coupling regime
> cα α the electron can spontaneously create from the va-
cuum the electron-hole pair (in the same way as the super-
critical charge creates electron-hole pairs). The initial elec-
tron attracts the hole and forms a bound state (an exciton)
and the emitted electron (which also has the supercritical
charge) can spontaneously create another pair, etc. The pro-
cess of creating pairs continues leading to the formation
of excitonic condensate and, as a result, the quasiparticles
acquire a gap. The semimetal-insulator transition in graphene
is similar to the chiral symmetry breaking phase transition
in strongly coupled QED studied in the 1970s and 1980s
(for a review see Ref. 34). The latter QED transition in-
duced by strong electromagnetic fields was searched in
experiments in heavy-ion collisions [35].
To stay closer to the experimental situation, one should
make a further step by considering electron states in the
field of two Coulomb centers, both like and unlike charg-
ed. The electron states in the field of charged impurities in
graphene and in the presence of a magnetic field are also of
considerable interest from the experimental point of view.
It was shown in Refs. 36, 37 that the strength of a charged
impurity can be tuned by controlling the occupation of
Landau-level states with a gate voltage.
All these topics are considered in the present review
paper which is organized as follows. The phenomenon of
the supercritical charge instability is briefly discussed in
Sec. 1. In Sec. 2, we analyse the electron states in gapless
and gapped monolayer graphene. The experimental data of
the observation of the atomic collapse in graphene are pro-
vided in Sec. 2.3. The impact of a magnetic field on the
supercritical charge problem in graphene is studied in
Sec. 3. Two-electron bound states of a Coulomb impurity
are considered in Sec. 4. The atomic collapse in the field of
two charged impurities is investigated in Sec. 5. The dipole
problem is studied in Sec. 6. The specifics of the atomic
collapse in bilayer graphene is considered in Sec. 7. The
results are summarized and conclusions are given in Sec. 8.
2. Atomic collapse in monolayer graphene
The electron quasiparticle states in the vicinity of the
K± points of graphene in the potential ( )V r of charged im-
purities are described by the following Dirac Hamiltonian
in 2 1+ dimensions:
( , ) = ( ),F zH Vξ + ξ∆σ +p p rv σ (1)
where Fv is the Fermi velocity of graphene, = i−p ∇ is the
canonical momentum, iσ are the Pauli matrices, ∆ is a
quasiparticle gap, and ξ is an index, which corresponds to
the valley K+ ( = 1ξ + ) or K− ( = 1ξ − ). Although the pris-
tine graphene is gapless, a quasiparticle gap ∆ can be gen-
erated if graphene sheet is placed on a substrate and two
carbon sublattices become inequivalent because of interac-
tion with the substrate (for band structure calculation of
such a configuration see, for instance, Ref. 38). The gap
can arise also in graphene ribbons due to geometrical quan-
tization [39] or due to many-body electron correlations
[24–33].
The Hamiltonian (1) acts on two component spinor sξΨ
which carries the valley ( = )ξ ± and spin ( =s ± ) indices.
We will use the standard convention: = ( , )T
s A B K s+ +
Ψ ψ ψ ,
whereas = ( , )T
s B A K s− −
Ψ ψ ψ , and ,A B refer to two sublat-
tices of hexagonal graphene lattice. Since the interaction
potential does not depend on spin, we will omit the spin in-
dex s in what follows. Further, for the sake of definiteness,
we will consider electrons in the K+ valley. The Hamilto-
nians at two valleys are related by means of the time rever-
sal operator 2 1= is KΘ σ :
1( , = 1) = ( , = 1),H H−Θ ξ + Θ − ξ −p p (2)
where 2s is the Pauli spin matrix and K is the complex con-
jugation.
The supercritical instability in the field of a single
charged impurity was studied quite in detail in the litera-
ture [12–15,21,40–44]. In this section we will summarize
its main features.
2.1. Resonance states in gapless graphene
in quasiclassical approach
Let us start our analysis with the case of gapless gra-
phene. Since massless particles cannot form bound states,
the atomic collapse is revealed for massless particles through
resonance states which appear when the Coulomb potential
strength exceeds a certain critical value = 1/ 2cα . In order
to demonstrate the presence of these states, it is instructive
to begin with the semiclassical analysis. We follow in this
subsection the derivation in Ref. 12.
In relativistic classical theory, the electron trajectories
can spiral around the charged center and eventually fall
down on it [45] if the electron angular momentum is small
enough 2< = e /cM M Z c .
These states can be constructed quasiclassically from
relativistic dynamics described by the Hamiltonian
= | | ( )FH V r+pv , where 2( ) = e / ( )V r Z r− κ and κ is a di-
electric constant. The collapsing trajectories with angular
momenta 2< = e / ( )c FM M Z κv are separated from non-
falling trajectories by a centrifugal barrier. This is mani-
fested in the expression for the radial momentum square
Low Temperature Physics/Fizika Nizkikh Temperatur, 2018, v. 44, No. 5 493
E.V. Gorbar, V.P. Gusynin, and O.O. Sobol
22 2
2 2
2
e= .r F
Z Mp E
r r
−
+ − κ
v (3)
Clearly, there is a classically forbidden region, the an-
nulus 1 2< <r r r , 2
1,2 = ( e / )/ | |Fr Z M Eκ v , where the
right-hand side of Eq. (3) is negative. The quasi-stationary
states trapped by this barrier are obtained from the Bohr–
Sommerfeld quantization
1
0
=
r
r
r
p dr nπ∫ ,
where 0r is a regularization parameter, which is of order of
lattice spacing. Evaluating the integral with logarithmic
accuracy, we obtain 2
0ln ( e / ( | |)) =Z r E nγ κ π , where
2 2 1/2( )cM Mγ ≡ − , which gives the quasi-Rydberg states
2
/
0
e e , > 0.n
n
ZE n
r
−π γ≈ −
κ
(4)
The energies of these states converge to zero, 0nE → ,
at large n, whereas their radii diverge, similar to the Ry-
dberg states in the hydrogen atoms. To find the transparen-
cy of the barrier, we integrate Im rp and obtain the tunnel-
ing action
( )
22 2
2
1
= = .
r
c
c
Fr
MM ES dr M
rr
− + π − γ
∫ v
(5)
Taken near the threshold 0γ ≈ , the transparency 2 /e S−
gives the width | | exp ( 2 )n nE ZΓ − π α , where 2= e /( )Fα κv
is the effective coupling constant. The quasi-Rydberg states
manifest themselves in the local density of states that can
be probed experimentally. Also, resonance scattering on
the quasi-bound states manifests itself in the dependence of
transport properties on the carrier density. For supercritical
potential strength | | > 1/ 2Zα there are oscillations of the
Ohmic conductivity which have a characteristic form of
Fano resonances centered at nE [12]. In this regime the
conductivity exhibits peaks at the densities for which the
Fermi energy FE equals nE . The peak position is highly
sensitive to the potential strength Zα, changing by an or-
der of magnitude when Zα varies from –1.0 to –1.3.
It is instructive to compare these results to the exact so-
lution of the Coulomb center problem that we do in the
next section.
2.2. Supercritical instability in graphene
2.2.1. Gapped graphene, subcritical regime
Now, let us include into consideration a quasiparticle
gap that on the one side makes more transparent the deri-
vation of the instability condition (diving of the lowest
energy level into the negative continuum), while on the
other hand takes into account a possible presence of a gap
due to the interaction with a substrate. In this subsection,
we follow the study performed in Ref. 21. The electron
quasiparticle states in graphene in the field of a single Cou-
lomb impurity are described by Dirac Hamiltonian (1) with
a regularized Coulomb potential
2 2
0 0
0
( ) = , ( > ), ( ) = , ( < )Ze ZeV r r r V r r r
r r
− −
κ κ
. (6)
As we discussed in the Introduction, to avoid the fall-to
center problem we should regularize the Coulomb potential
at small distances. Potential (6) represents the simplest
“cutoff” regularization. Since the Hamiltonian (1) with
potential (6) commutes with the total angular momentum
operator
= =
2z z z zJ L S i ∂
+ − + σ
∂φ
,
we seek eigenfunctions in the following form:
( 1/2)
( 1/2)
e ( )1= .
e ( )
i j
i j
a r
r i b r
φ −
φ +
Ψ
(7)
Then we obtain a system of two coupled ordinary dif-
ferential equations of the first order
( )( 1/ 2) = 0,
( )( 1/ 2) = 0.
F
F
a E V ra j b
r
b E V rb j a
r
+ ∆ −′ − + +
− ∆ −′ + − −
v
v
(8)
It is convenient to define the quantities = / ,FEε v
= / Fm ∆ v , and 2= / = e /g Fα α κ κv .
The discrete spectrum of Eqs. (8) exists for | | < mε . In
this case it is convenient to define
2 2= , = 2 , = ( ),
2
= ( )
2
mu m ur a g f
mb g f
+ ε
− ε ρ −
− ε
+
(9)
and rewrite Eqs. (8) in the region r > r0 as follows:
1 = 0,
2 2
1 = 0.
2 2
mg g Z f j Z
u u
mf f Z g j Z
u u
ρ ε ′ρ + − − α + + α
ρ ε ′ρ − + − α + − α
(10)
Substituting f from the first equation into the second
one, we obtain the equation for the g component
2 2 2
2
2 2
1 1
1 2 4 = 0,
4
Z j Zd g u g
d
ε + α − + α
+ − + +
ρρ ρ
(11)
494 Low Temperature Physics/Fizika Nizkikh Temperatur, 2018, v. 44, No. 5
Electron states in the field of charged impurities in two-dimensional Dirac systems
which is the well-known Whittaker equation [46]. Its gene-
ral solution is
1 , 2 ,= ( ) ( ),g C W C Mµ ν µ νρ + ρ 1= ,
2
Z
u
αε
µ +
2 2 2= .j Zν − α
(12)
Taking into account the asymptotic of the Whittaker func-
tions , ,( ), ( )k kW z M zν ν at infinity,
, ( ) e (2 ) ,urW ur− µ
µ ν ρ
,
(1 )( ) e (2 ) ,
1
2
urM ur −µ
µ ν
Γ + ν
ρ
Γ −µ + ν
,r →∞ (13)
we find that the regularity condition at infinity requires
2 = 0C . Then the first equation in (10) gives the following
solution for the f component in the region II ( 0>r r ):
1 1 ,
2
= ( ).II Z
u
mf C j Z W
u ε
− + α ν
− α ρ
(14)
Solutions in the region I ( 0<r r ) are easily obtained
2
2
1 | 1/2|
0
= ,I j
Zb A rJ r m
r+
α ε + −
(15)
2
20
1 | 1/2|
0 0
/
= sgn ( ) ,
/I j
Z r m Za A j rJ r m
Z r m r−
ε + α + α ε + − ε + α −
(16)
where 1A is a constant and we took into account the infra-
red boundary condition which selects only regular solution
for Ib and Ia . Energy levels are determined through the
continuity condition of the wave function at 0=r r ,
= =0 0
= ,I II
I IIr r r r
b b
a a
(17)
that gives the equation
1 ,
2
= 0
1 ,
2
( )
1| = ,
1( )
Z
u
r r
Z
u
W
k
Z m kj W
u
αε
+ ν
αε
− + ν
ρ
+
α − − ρ
| 1/2|0
0 | 1/2|
( )/
= sgn ( ) ,
/ ( )
j
j
JZ r mmk j
u Z r m J
+
−
ρε + α −+ ε
ε + α + ρ
(18)
2 2 2
0 0= ( ) .Z r m rρ α + ε −
We analyze this equation in the limit 0 0r → where we
can use the asymptotical behavior of the Whittaker func-
tion at 0ρ→ ,
1 1
2 2,
(2 ) ( 2 )( ) .
1 1( ) ( )
2 2
W
−ν +ν
µ ν
Γ ν Γ − ν
ρ ρ + ρ
Γ −µ + ν Γ −µ −ν
(19)
In the limit 0 0r → Eq. (18) reduces to the following one,
2
0
1
( 2 ) (2 ) =
(2 ) 1
Z
u ur
Z
u
ν
ε Γ + ν − α Γ − ν
εΓ ν Γ − ν − α
0
0
0
( ) ( )
= ( ),
( ) ( )
Z m Z mj k j
u u O r
Z m Z mj k j
u u
α + ε α − ε + ν − + − ν −
− +
α + ε α − ε − ν − + + ν −
(20)
where
| 1/2|
0
| 1/2|
( )
= sgn ( ) ( , ).
( )
j
j
J Zm mk j Z j
u J Z u
+
−
α+ ε + ε
≡ σ α
α
(21)
Equation (20) can be rewritten in more convenient form
2
0
1
( 2 ) (2 ) =
(2 ) 1
Z
u ur
Z
u
ν
ε Γ + ν − α Γ − ν
εΓ ν Γ − ν − α
( )
( , )= .
( ) ( , )
Z mj j Z Z ju
Z m j Z Z jj
u
α − ε
− ν − + ν − ασ α
−
α − ε − ν − ασ α+ ν −
(22)
In the limit 0 0r → the energy levels are determined by the
poles of the gamma function (1 / )Z uΓ + ν − αε and by a zero
of the right hand side of Eq. (22), this leads to the familiar
result (analogue of the Balmer’s formula in QED) [47] (re-
derived also in [40]),
1/22 2
, 2
= 0,1,2,3,..., > 0,
= 1 ,
= 1,2,3,..., < 0.( )
n j
n jZm
n jn
−
α
ε +
ν +
(23)
The bound states for 1n ≥ are doubly degenerate,
, ,=n j n j−ε ε . The lowest energy level is given by
2
0, =1/2 = 1 (2 ) .j m Zε − α (24)
If Zα exceeds 1/ 2, then the ground state energy (24)
becomes purely imaginary, i.e., the fall-to-center phenome-
non occurs [12,13,43,44]. In fact, all energies ,1/2nε be-
come complex for > 1/ 2Zα . The unphysical complex
energies indicate that the Hamiltonian of the system is not
a self-adjoint operator for supercritical values > 1/ 2Zα
and should be extended to become a self-adjoint operator.
According to [2,3], nonzero 0r resolves this problem. For
>Z jα , ν is imaginary for certain j and for such j we
Low Temperature Physics/Fizika Nizkikh Temperatur, 2018, v. 44, No. 5 495
E.V. Gorbar, V.P. Gusynin, and O.O. Sobol
denote = ,iν β 2 2 2= Z jβ α − . For finite 0r discrete levels
also exist for > 1/ 2Zα . Their energy decreases with in-
creasing of Zα until they reach the lower continuum. The
behavior of lowest energy levels with = 1/ 2j as functions
of the coupling Zα is shown in Fig. 1(a).
The critical charge cZ that corresponds to diving into
the continuum is obtained from Eq. (22) setting = mε −
there and using the corollary of the Stirling formula:
2 ln( ) e ,
( )
iy xx iy x
x iy
Γ +
→ → +∞
Γ −
.
We come at the equation
2 ln (2 )0 ( , ) (1 2 )e = ,
( , ) (1 2 )
i Z mr i j Z Z j i
i j Z Z j i
− β α β − + ασ α Γ − β
− β− + ασ α Γ + β
(25)
or,
0ln (2 ) = arg ( ( , ) )Z mr Z Z j j i−β α ασ α − + β +
arg (1 2 ) , 0, 1, ...i n n+ Γ − β + π = (26)
It is not difficult to check that for = 1/ 2j and = 1n the
critical coupling cZ α approaches the value 1/2 for 0 0mr → .
The dependence of the critical coupling cZ α on 0mr for
= 1/ 2j is shown in Fig. 1(b).
2.2.2. Gapped graphene, supercritical regime
Let us analyze Eq. (18) in the supercritical case
> 1/ 2Zα and show that there are resonant states for
| | > mε (we define the gap > 0∆ ). The Whittaker function
, ( )Wµ ν ρ with = 1/ 2 /Z uµ + αε , 2 2 2= j Zν − α describes
bound states for | |< mε which are situated on the first
physical sheet of the variable u and for which Re > 0u
(see, Eq. (14)). The quasistationary states are described by
the same function , ( )Wµ ν ρ and are on the second unphysi-
cal sheet with Re < 0u . We shall look for the solutions
corresponding to the quasistationary states which define
outgoing hole waves at r →∞ with
Re < 0, Im < 0, Re < 0, Im < 0.u uε ε (27)
For solutions with 2 2 2>Z jα resonance states are deter-
mined by Eq. (18) for bound states where ν is replaced by
= iν β . We will consider the states with = 1/ 2j which
correspond to the 1/2nS -states, in particular, the lowest en-
ergy state belongs to them. The corresponding equation
then takes the form
1 ,
2
= 0
1 ,
2
0 1
0 0
2 2 2
0 0
( )
1| = ,
1 1( )
2
/ ( )
= ,
/ ( )
= ( ) .
Z i
u
r r
Z i
u
W
k
Z m kW
u
Z r m Jmk
u Z r m J
Z r m r
αε
+ β
αε
− + β
ρ
+
α − − ρ
ε + α − ρ+ ε
ε + α + ρ
ρ α + ε −
(28)
The analytical results can be obtained for the near-
critical values of Z when 1/ 2 1Zα − . We assume that
0| 2 | 1ur , then using the asymptotic of the Whittaker
function, we find
2
0
1
(1 2 )(2 ) =
(1 2 ) 1
i
Zi
i uur
Zi i
u
β
αε Γ + β− Γ − β
αεΓ + β Γ − β−
1
0
1
0
( )11 ( )
2 ( )2= .
1 ( ) ( )1
2 2 ( )
J ZZ m i Zi J Zu
Z m J Zi i Z
u J Z
αα − ε + β− α− β− α
α − ε α
+ β− − β− α
α
(29)
Expanding Eq. (29) in the near critical region in powers of
2 2= 1/ 4Zβ α − , we find the following equation:
Fig. 1. (Color online) The lowest energy levels as functions of
Zα . Red lines correspond to the pure Coulomb potential (they
exist only for < 1 / 2Zα ); black solid lines are numerical solution
for = 1 / 2,j 0 = 0.01mr ; black dashed line is a numerical solu-
tions for = 1 / 2,j − 0 = 0.01mr (a). The critical coupling as
a function of 0mr for the 1/21S level (b).
496 Low Temperature Physics/Fizika Nizkikh Temperatur, 2018, v. 44, No. 5
Electron states in the field of charged impurities in two-dimensional Dirac systems
2 2 2 0
0
0 1
(1 / 2)( 2 ) = 1 4
(1 / 2) (1 / 2)
i Ji m r i
J J
β
− ε − + + β + −
2 2
1 1(1) 1
2 2
1
i
mm i
m
ε + Ψ − Ψ − −
ε −ε − + ε +
. (30)
Here ( )xΨ is the psi-function and we put 2 2=u i m− ε −
where 2 2Im < 0mε − on the second sheet.
It is instructive to consider resonant states in the vicini-
ty of the level = mε − when bound states dive into the low-
er continuum and determine their real and imaginary parts
of energy. First of all, nonzero m increases the value of
the critical charge. Indeed, using Eq. (26), we obtain that
the critical value cZ α for = 1/ 2j scales with m like (see
Fig. 1(b))
2
2
0
0
0 1
1 ,
2 ( )ln
2 (1/ 2)
= exp 2 (1) 0.21.
(1/ 2) (1/ 2)
cZ
cmr
J
c
J J
π
α +
− Ψ − ≈ −
(31)
Note that the dependence of the critical coupling on 0mr is
quite similar to that in the strongly coupled QED [34,48].
For > cZ Z , using Eq. (30), we find the following reso-
nant states:
/ 23 3= 1 e , = ,
8 8
b c
c
m b i b−π β −βπ π ε − + + ββ
(32)
where 2= ( ) 1/ 4c cZβ α − . Like in QED [49] the imagi-
nary part of energy of these resonant states vanishes expo-
nentially as cZ Z→ . Such a behavior is connected with tun-
neling through the Coulomb barrier in the problem under
consideration. For the quasielectron in graphene in a central
potential ( )V r , expressing the lower component of the Di-
rac spinor (7) through the upper one and following [3,49],
we obtain an effective second order differential equation in
the form of the Schrödinger equation
2( ) ( ) ( ) = 0,'' r k r rχ + χ
1 1( ) = exp ( ).
2
'Va r dr r
r m V
− χ ε + −
∫
(33)
Here
2 2
2 ( ) = 2( ( )), = , = ,
2 F
m Vk r U r Vε −
−
v
(34)
and we represent the effective potential as the sum of two
terms 1 2=U U U+ , where 1U is the effective potential for
the Klein–Gordon equation and 2U takes into account the
spin dependent effects,
2
1 2
( 1)= ,
2 2
V j jU V
r
−
ε − +
(35)
2
2
1 3 2= .
4 2 ( )
'' ' 'V V jVU
m V m V r m V
+ + ε + − ε + − ε + −
(36)
Note that Eq. (33) and the potentials (35), (36) coincide
with the corresponding equations in QED [3]. One can
show that in the near-critical regime ( cZ Z→ , = 1/ 2j ,
and = mε − ) the effective potential ( )U r has the Coulomb
barrier (see Fig. 3 below), which prevents the delocaliza-
tion of the wave function.
The tight-binding approach (solved exactly by using nu-
merical techniques) was compared with the continuum ap-
proach based on the Dirac equation in Ref. 14. It was shown
that the latter provides a good qualitative description of
the problem at low energies when properly regularized. On
the other hand, the Dirac description fails at moderate to
high energies and at short distances when the lattice de-
scription should be used.
2.2.3. Gapless graphene
We consider now the case of gapless graphene, = 0m .
Writing = | | eiγε ε Eq. (30) takes the form
0ln (2 | | )
2
r i π ε + γ −
0
0 1
(1/ 2) 1 12 (1) 1 ,
(1/ 2) (1/ 2) 2 2 1
J i n
J J i
π + Ψ − Ψ − − − − + β
= 1, 2, .n (37)
We find
(0) 1
0 2 2
= e exp =
1/ 4
i
n
nar
Z
− γ π
ε −
α −
1
0 2 2
= (1.18 0.17 ) exp , = 1, 2, ,
1/ 4
ni r n
Z
− π
− + −
α −
(38)
where
= 1 coth 3.28,
2 2
π π γ + ≈
0
0 1
2 (1/ 2)1= exp 2 (1) 1 Re 1 1.19.
2 (1/ 2) (1/ 2) 2
J ia
J J
+ Ψ − − Ψ − ≈ −
(39)
These results are in agreement with Eq. (4) and Refs. 12, 13.
The energy of quasistationary states (38) has a characteris-
tic essential-singularity type dependence on the coupling
constant reflecting the scale invariance of the Coulomb po-
tential. The infinite number of quasistationary levels is re-
lated to the long-range character of the Coulomb potential.
Note that a similar dependence takes place in the supercri-
tical Coulomb center problem in QED [50].
Since the “fine structure constant” 2 / 2.2Fe ≈v in gra-
phene, an instability could potentially appear already for
the charge = 1Z . However, in the analysis above we did
Low Temperature Physics/Fizika Nizkikh Temperatur, 2018, v. 44, No. 5 497
E.V. Gorbar, V.P. Gusynin, and O.O. Sobol
not take into account the vacuum polarization effects. Con-
sidering these effects and treating the electron-electron
interaction in the Hartree approximation, it was shown in
Ref. 42 that the effective charge of impurity effZ is such
that the impurity with bare charge = 1Z remains subcriti-
cal, 2
eff e /( ) < 1/ 2FZ κv , for any coupling 2 /( )Fe κv ,
while impurities with higher Z may become supercritical.
For finite m and in the case | | ,mε Re < 0ε , expand-
ing Eq. (30) in /m ε we get up to the terms of order 2 2/m ε ,
2 2
(0)
2= 1 (0.29 0.23 ) ,
2 n
m m m i
ε − ε − + − ε ε ε
= 1, 2,n . (40)
The resonant states with (0)
nε describe the spontaneous
emission of positively charged holes when electron bound
states dive into the lower continuum in the case = 0m . In
order to find corrections to these energy levels due to non-
zero m, we seek solution of Eq. (40) as a series ( )
=0
= k
k
∞
ε ε∑
with ( )kε of order km and easily find the first two terms
2
(0)
(0)= (0.24 0.20 ).
| |
n n
n
mm iε ε − + +
ε
(41)
Since (0)Im < 0nε , the appearance of a gap results in de-
creasing the width | Im |ε of quasistationary states and, there-
fore, increases stability of the system. Also, as we showed
above, the critical value 0( )cZ mrα , determined by the con-
dition of appearance of a nonzero imaginary part of the
energy, increases with the increase of m. Thus there are
two possibilities for the system with supercritical charge to
become stable: to create spontaneously electron-holes pairs
and shield the charge or to generate spontaneously the
quasiparticle gap. In the problem of the supercritical Cou-
lomb center only the first possibility can be realized, which
is already due to the formulation of the problem as the one-
particle one. The second possibility — dynamical genera-
tion of the gap — was studied in Refs. 21–33.
Considering the many-body problem of strongly inter-
acting gapless quasiparticles in graphene, it was shown that
the Bethe-Salpeter equation for an electron-hole bound
state contains a tachyon in its spectrum in the supercritical
regime > cα α , the critical constant = 1.62cα in the static
random-phase approximation [21] and = 0.92cα in the
case of the frequency-dependent polarization function [28].
The tachyon states play the role of quasistationary states in
the problem of the supercritical Coulomb center and lead
to the rearrangement of the ground state and the formation
of excitonic condensate. Thus, there is a close relation be-
tween the two instabilities, in fact, the tachyon instability
can be viewed as the field theory analog of the fall into the
center phenomenon and the critical coupling cα is an ana-
log of the critical coupling cZ α in the problem of the Cou-
lomb center. The physics of two instabilities is related to
strong Coulomb interaction.
2.3. Experimental observation of the atomic collapse
in graphene
Univalent charge impurities, such as K, Na, or 3NH , all
commonly used in graphene, are on the border of the su-
percritical regime. To investigate this regime experimental-
ly, one can use divalent or trivalent dopants such as alka-
line-earth or rare-earth metals. However, the observation of
atomic collapse in the field of supercritical impurities has
remained elusive for some time due to the difficulty of
producing highly charged impurities.
For the first time, the supercritical Coulomb behavior
was observed in atomically-fabricated “artificial nuclei”
assembled on the surface of a gated graphene device in
Ref. 16. Calcium atoms were deposited onto the graphene
device at low temperature < 10 KT . Then graphene was
warmed up before returning to lower temperature, thus
causing the Ca adatoms to thermally diffuse and bind into
dimers. Further, as charges are transferred from a Ca dimer
into graphene band states, the Ca dimer becomes positively
charged. By making use of the density functional theory
calculations, it was found that Ca dimers acquire an effec-
tive positive charge 0.4e.
The tunable charge centers were synthesized by pushing
together Ca dimers using the tip of a scanning tunneling
microscope (STM) (see insets to Fig. 2(a)–(c)), thus allow-
ing creation of supercritical Coulomb potentials from sub-
critical charge elements. The scanning tunneling spectros-
copy was used to observe the emergence of atomic-collapse
electronic states extending further than 10 nm from the cen-
ter of artificial nuclei in the supercritical regime ( > cZ Z ).
Here, the effective charge Z is defined as the screened
cluster charge where the effects of intrinsic screening due
to graphene band polarization and the substrate are taken
into account, and the critical value is 2= / 2e 0.25c FZ v .
By tuning the graphene Fermi level FE via electrostatic
gating the atomic collapse behavior was observed.
Experimentally, the local density of states (LDOS) is
measured by means of the STM technique. A sharp STM
tip scans over a graphene piece and measures the electric
current I from the surface due to the tunneling effect. This
current depends on the voltage V between tip and a sample
and its derivative with respect to V is proportional to
the LDOS, / ( , )dI dV D E r where ( , )D E r is given by
Eq. (64) below. The curves in Fig. 2 show the differential
conductance /dI dV (and thus the LDOS) as a function of
the bias voltage V , hence the energy =E eV . The various
curves in panels (a–c) correspond to different distances
from the charge center in the range of about 2–20 nm.
Spectra acquired near 1-dimer clusters (Fig. 2(a)) dis-
played electron-hole asymmetry as well as an extra oscilla-
tion in the LDOS at high energies above the Dirac point.
For the 4-dimer cluster, the resonance is clearly observed
498 Low Temperature Physics/Fizika Nizkikh Temperatur, 2018, v. 44, No. 5
Electron states in the field of charged impurities in two-dimensional Dirac systems
close to the Dirac point (Fig. 2(b)). For the 5-dimer cluster,
the resonance shifts below the Dirac point (Fig. 2(c)). The
formation of this resonance (or quasi-bound state) as nu-
clear charge increases is the “smoking gun” for the atomic
collapse. The experimental data suggest that clusters with
just one or two Ca dimers are in the subcritical regime. The
clusters composed of four or more dimers are either (for
four dimers) transitioning into or (for five dimers) have
fully entered the supercritical regime, as evidenced from
panels (b) and (c) in Fig. 2. For these clusters, / cZ Z is
determined by matching the quasi-bound state resonance
energy between the simulation and experiment. The main
features seen in the experimental data are well reproduced
by the Dirac equation simulations in Ref. 16.
In order to check that the magnitude of / cZ Z extracted
for Ca dimers from the Dirac equation fits is physically
reasonable, a completely separate density functional theory
calculation of the charge state expected for a Ca dimer
adsorbed to graphene was performed [16]. This calculation
(which had no fitting parameters) yielded a single-dimer
charge ratio of / = 0.6 0.3cZ Z ± . This is in agreement
with the value / = 0.5 0.1cZ Z ± obtained via Dirac equa-
tion simulations, and thus lends further support to overall
interpretation of the data. The behavior of the quasi-bound
state observed for high-Z artificial nuclei depends on
whether it is occupied by electrons or empty. For the de-
tails of this doping dependence see the original paper [16].
3. Supercritical instability in a magnetic field
As we discussed in the Introduction, the supercritical
charge instability in a many-body system leads to much
more dramatic consequences compared to the single-particle
problem of the Coulomb center. Like the Cooper instability
in the theory of superconductivity, the QED supercritical
coupling instability is resolved only through the formation
of a condensate of electron-positron pairs generating a mass
gap in the spectrum [34]. It was shown in [54–57] that
magnetic field catalyses gap generation in relativistic-like
systems and even the weakest attraction leads to the for-
mation of a symmetry breaking condensate. Therefore, the
many-body system is always in the supercritical regime
once there is an attractive interaction. The magnetic cataly-
sis plays an important role in quantum Hall effect studies
in graphene [58], where it is responsible for lifting the de-
generacy of the Landau levels.
In QED in (3+1) dimensions, the Coulomb center prob-
lem in a magnetic field was studied for massive fermions
in [59,60]. There it was found that the magnetic field con-
fines the transverse electronic motion and the electron in a
magnetic field is closer to the nucleus than in the case
where magnetic field is absent. Thus, it feels stronger Cou-
lomb field. Therefore, cZ α decreases with B. The Dirac
equation for (2+1)-dimensional quasiparticles in graphene
in the Coulomb potential in a magnetic field was consid-
ered in Ref. 61 where exact solutions were found for cer-
tain values of magnetic field, i.e., this problem furnishes an
example of the so-called quasi-exactly solvable models.
However, no instability or resonance was found.
We would like to stress that the presence of a constant
magnetic field changes qualitatively the supercritical Cou-
lomb center problem. Indeed, if magnetic field is absent,
then the supercritical Coulomb center instability leads to a
resonance which describes an outgoing positron propagat-
ing freely to infinity. However, since charged particles in
a plane perpendicular to a magnetic field do not propagate
freely to infinity, such a behavior is impossible for the in-
plane Coulomb center problem in graphene in an out-of-
plane magnetic field. Therefore, a priori it is not clear how
the supercritical instability manifests itself in the Coulomb
center problem in a magnetic field. This question was stud-
ied in Ref. 62. We would like to note that the role of
a magnetic field for the atomic collapse in graphene is ra-
ther subtle and different conclusions on this issue were
drawn in the literature [63–65].
Fig. 2. (Color online) Evolution of charged impurity clusters from subcritical to supercritical regime. (a)–(c) /dI dV spectra measured at
different distances from the center of Ca-dimer clusters (i.e., artificial nuclei) composed of 1, 4, and 5 dimers. “Center” here is defined
as the average coordinate of dimers within a cluster. All spectra were acquired at the same back-gate voltage ( = 30 VgV − ) and each was
normalized by a different constant factor to account for exponential changes in conductivity due to location-dependent tip-height changes
[51–53]. Insets: STM topographs of atomically fabricated Ca-dimer clusters. The nuclear charges / cZ Z of (a) 0.5 (1-dimer cluster),
(b) 1.8 (4-dimer cluster), and (c) 2.2 (5-dimer cluster). Black dashed lines indicate Dirac point, red arrows indicate atomic collapse state
observed in experiment. This figure is an adapted version of the corresponding figure from Ref. 16.
Low Temperature Physics/Fizika Nizkikh Temperatur, 2018, v. 44, No. 5 499
E.V. Gorbar, V.P. Gusynin, and O.O. Sobol
In the presence of a charged impurity, degenerate Lan-
dau levels convert into bandlike structures due to lifting the
orbital degeneracy. For zero chemical potential, as the
charge of impurity increases, the energy level with the
quantum numbers = 0n , = 1/ 2j − comes close to the
highest energy state of the level = 1n − . In the absence of
magnetic field, the corresponding bound state would dive
into the lower continuum and further increase of the charge
of impurity would produce a resonance. The situation is
qualitatively different in the presence of a magnetic field as
the energy curves with the same momenta j never cross.
The results clearly demonstrate this phenomenon of the
level repulsion between the sublevels with the same j and
the formation of a quasiresonance state when the impurity
charge exceeds a critical value. In such a case we observe a
redistribution of profiles of radial distribution functions
with the same orbital momentum among lower Landau
levels 1n ≤ − .
3.1. The Coulomb center in a magnetic field
Let us consider the electron states in gapped graphene
with a single charged impurity in a magnetic field. The cor-
responding Hamiltonian could be obtained from Eq. (1) by
the standard substitution = e
ci→ − ∇+p p A , where
< 0e− is the electron charge and the vector potential
= / 2( , )B y x−A in the symmetric gauge describes magnet-
ic field perpendicular to the plane of graphene. We regular-
ize the Coulomb potential of an impurity by introducing
a parameter 0r of the order of the graphene lattice spacing.
Then the regularized interaction potential of the impurity
with charge Ze is given by
2
2 2
0
e( ) = ZV
r r
−
κ +
r . (42)
It is convenient to introduce the magnetic length
= / | |Bl c eB and the dimensionless quantity 2= e /( )FZζ κv
which characterizes the strength of the bare impurity. Since
the total angular momentum is conserved, we use the polar
coordinates ( , )r φ and seek eigenfunctions in the form (7).
Then the Dirac equation takes the form
2
2
1/ 2 ( ) = 0,
2
1/ 2 ( ) = 0.
2
FB
FB
j r E V ra a a b
r l
j r E V rb b b a
r l
+ + ξ∆ − ′ − − +
− − ξ∆ − ′ + + −
v
v
(43)
Eliminating, for example the function ( )a r , one can get
a second order differential equation for ( )rχ defined by the
relation
1/2 ( )[ ( )] ( ) = .b rE V r r
r
− ξ∆ − χ (44)
We obtain the following Schrödinger-like equation:
( ) ( ) ( ) = ( ),r U r r r′′−χ + χ χ (45)
where
2 2
2( )F
E − ∆
=
v
(46)
and the effective potential, 1 2=U U U+ , reads
2
1 2 2 4 2
(2 ) ( 1) 1/ 2= ,
( ) 4F B B
V E V j j r jU
r l l
− + −
+ + +
v
(47)
2
2 2
1 3 2= .
2 2 ( )2 B
V V j r VU
E V E V r E Vl
′′ ′ ′ + − + − ξ∆ − − ξ∆ − − ξ∆ −
(48)
We plot the effective potential ( )U r near the K− point
( = 1ξ − ) for =E −∆ and = 1/ 2j − in Fig. 3, where the
energy barrier in the absence of magnetic field is clearly
seen, which leads to the appearance of resonances for suf-
ficiently large charge. We note that the equations for spinor
components ( )a r and ( )b r at the K− point can be obtained
from the equations in Sec. 2.2.1 at the K+ point by inter-
changing a b↔ and changing j j→ − since two points are
related by means of the time reversal transformation,
=K K− +
Ψ ΘΨ , introduced in Sec. 2. The presence of a mag-
netic field changes the asymptotic of the effective potential
at infinity and, thus, forbids the occurrence of resonance
states. This feature distinguishes qualitatively the Coulomb
center problem in a magnetic field from that at = 0B .
Unfortunately, Eq. (45) belongs to the class of equa-
tions with two regular and one irregular singular (at =r ∞)
points, and cannot be solved in terms of known special
functions. In the regime 0Zα → , we can find it using per-
turbation theory. For = 0Zα , the corresponding solutions
are the well known Landau states degenerate in the total
angular momentum j . The Coulomb potential of impurity
removes degeneracy in j and the eigenenergies split into
series of sublevels resulting in an j dependent energy njE .
The energy downshift is largest for 0 jE and diminishes
with increasing –j. For the = 0n level with (0) (0)
0 0= =jE E ∆
the normalized wave function has the form (at the K− point)
Fig. 3. The potential ( )U r as a function of a distance from the
Coulomb center at zero and nonzero magnetic field, near the K−
point for =E −∆ and = 1 / 2j − .
500 Low Temperature Physics/Fizika Nizkikh Temperatur, 2018, v. 44, No. 5
Electron states in the field of charged impurities in two-dimensional Dirac systems
2 2 /2/4 2
2
0
( 1)( , ) = e ,
2 ! e
2
M Mr lBM iM
B
B
r rl M
l
−
− φ
−
Ψ φ
π
(49)
where = ( 1/ 2) = 0, 1, 2, ...M j− + is the orbital quantum
number. Energy corrections Mε of perturbed states of the
Landau level (0)
0 =E ∆ are found from the secular equation
1 2
| | = 0M MVε − , (50)
where 1 2M MV is a matrix element of the potential on states
(49). Since
1 2M MV is a diagonal matrix, we easily obtain
22 2 1 /2
2 2
0 0
= = =
!2
M
M MM M
B
Ze d eV
M l
∞ + −ρρρ
ε −
κ ρ +ρ
∫
( )
2
2 1 2
0 01= 1 ,3 / 2 ; / 2 ,
2
M
M
B
Ze M M
l
+
+
− ρ Ψ + + ρ
κ
(51)
where ( , ; )a c zΨ is the confluent hypergeometric function,
0 0= / Br lρ . For small 1Zα we can use the unregularized
Coulomb potential, then setting 0 = 0ρ in Eq. (51) we get
2 1( )
2= .
2 ( 1)M
B
Ze M
l M
Γ +
ε −
κ Γ +
(52)
Thus at large M the energy levels accumulate near the value
=E ∆ ,
2
(0)
0 0= .
2M M
B
ZeE E
l M
+ ε ∆ −
κ
(53)
The largest correction by modulus 0 = / 2F BZ lε − α πv
is for the state with = 0M . Naturally, in the lowest order
of perturbation theory, the energy linearly decreases with
the increase of the impurity charge. The numerical solution
of Eq. (45) shows that this behavior changes when the
charge exceeds a certain critical value and after that the
level repulsion occurs (see Fig. 4(a)).
For finite ∆ one can define the critical charge by the
condition (0)
00= =E E + ε −∆ when the lowest energy emp-
ty level descending from the upper continuum crosses the
energy level of a filled state. In the regime of small cou-
pling, 1Zα and 1Bl∆ , this gives
2 2
= .B
c
F
lZ ∆
α
πv
(54)
Clearly, this critical charge tends to zero as 0∆ → , while
the state with = 0M of the zero Landau level moves be-
low zero energy for any small impurity charge (its energy
is 0 = / 2F BZ lε − α πv ). The states connected with the
zero Landau level play an important role in the many-body
problem, e.g., in the formation of the excitonic condensate
and gap generation for quasiparticles [25,66] due to the
magnetic catalysis. In the case of a charged impurity in a
magnetic field, the negative energy states are filled and it is
physically more sensible to connect the critical charge with
the anticrossing of Landau levels in the negative energy re-
gion (see the discussion below).
Although, in view of the magnetic catalysis [56], a non-
zero gap is always generated in graphene in a perpendi-
cular magnetic field [24–27], this gap is rather small for
realistic magnetic fields. Therefore, it makes sense to ne-
glect it and see how levels with the same j evolve. Let us
solve Eqs. (43) numerically by using the shooting method.
In order to utilize this method, one should determine
the appropriate asymptote of the solution at 0r → for
2
0| ( ) | | (0) | = /( ) | |V r V Ze r E≈ κ
. At the K+ point it is
convenient to introduce the orbital quantum number
= 1/ 2 = 1m j M− − − . Then, for 0m ≥ ( 1/ 2)j ≥ the upper
component of (7) dominates and the leading behavior is
1( ) ma r r +
, for < 0m ( < 1/ 2j ) the lower component
dominates with ( ) mb r r− (see Ref. 67).
The numerical integration of Eq. (43) proceeds as fol-
lows. We take a “shot” from = 0r at a fixed value of ener-
gy solving the differential equations with the correct initial
Fig. 4. (Color online) The colormap of the LDOS at the impurity position ( = 0r ) as a function of coupling ζ and energy E in the mag-
netic field = 10B T. Black labels indicate the Landau level numbers n and orbital quantum numbers m (a). Critical coupling constant as
a function of magnetic field for two types of regularized impurity potential: displaced impurity, Eq. (42) — blue lines; cut-off potential
(6) — red lines. Solid lines correspond to anticrossing of the ˆ =H EΨ Ψ , = 0m and = 1n − , = 0m levels, dashed lines correspond to
anticrossing of the = 1n , = 1m − and = 1n − , = 1m − levels. In both panels the regularization parameter is chosen as 0 = 0.5r nm (b).
Low Temperature Physics/Fizika Nizkikh Temperatur, 2018, v. 44, No. 5 501
E.V. Gorbar, V.P. Gusynin, and O.O. Sobol
conditions and check the behavior of the wave functions at
r →∞. The latter may tend to +∞ for some values of energy
or to −∞ for other values. A physical solution is the solu-
tion for which the exponentially growing behavior of the
absolute value is absent. We find the corresponding value
of the energy of this solution by using the method of bisec-
tions. In all numerical calculations, we use 0 = 0.05 Br l .
The magnetic field modifies the energy spectrum of
electrons in the Coulomb field of the charged impurity
making all continuum states discrete and provides an effec-
tive scale given by the magnetic length. On the other hand,
the charged impurity removes the orbital degeneracy of
Landau levels transforming the latter into bandlike struc-
tures. Figure 4(a) shows the colormap of the LDOS at the
impurity position as a function of coupling ζ and dimen-
sionless energy /B FEl v in the magnetic field = 10B T.
Red lines correspond to Landau levels split into sublevels
with different orbital numbers. At the beginning, the curves
decrease linearly in accordance with Eq. (52). As the charge
of impurity increases, the curves, which correspond to the
= 0,1n Landau levels come close to lower curves, which
form a “quasicontinuum”. In the absence of magnetic field,
with further increase of the charge of impurity the corre-
sponding bound state would dive into the lower continuum
producing a resonance.
According to Fig. 4(a), the situation is qualitatively dif-
ferent in the case of zero gap when a magnetic field is pre-
sent as the curves with the same orbital number m never
cross each other. Instead, typical level repulsions are real-
ized (the well-known avoided crossing theorem [68] for-
bids a level crossing for two states with the same sym-
metry). We clearly see the repulsion between the levels
= 1n , = 1m − and = 1n − , = 1m − , as well as between the
levels = 0,n = 0m and = 1n − , = 0m . States with differ-
ent quantum numbers m simply cross each other without
repulsion. The situation is similar to that of a quantum
electrodynamical system of finite size [4,69]. In Fig. 4(b)
we plot the dependence of the critical charge on a magnetic
field B defined as the anticrossing points of the Landau
levels = 0n and = 1n − with = 0m , as well as of the levels
= 1n and = 1n − with = 1m − , in the negative energy re-
gion. The corresponding dependence at zero gap = 0∆ can
be very accurately fit by the following function:
cr 2
0
1= ,
2 ln ( / )B
A
cr l
ζ + (55)
where A and c are fitting parameters, which can be deter-
mined numerically. For the displacement regularization,
these parameters are 1 = 15.75A , 1 = 0.305c and, for the
cut-off regularization, they are 2 = 13.94A , 2 = 0.25c . The
increasing magnetic field strength causes the anticrossings
to appear at higher charge ζ in accordance with the obser-
vation in Ref. 65. The dependence of the critical charge
crζ on a magnetic field is similar to its dependence on a
gap in the absence of the field (see Eq. (31)).
Figure 5 shows the radial distribution function
2( ) = 2 | |nmW r rπ Ψ for the = 0m and = 0, 1, 2n − − states
for the three values of the impurity charge = 0.7ζ , 1.3, and
1.9. The second value corresponds to the states in the vi-
cinity of the avoided crossing. For a small charge of the
impurity (a), the electron density is weakly affected by the
impurity and the radial distribution functions of the above
mentioned states have one, two, and three maxima, respec-
tively. As the impurity charge increases, all leftmost max-
ima in ( )W r move to the impurity position = 0r and attain
their maximal values at 1.3ζ ≈ (b). In addition, a new max-
imum appears on the blue solid curve (as well as additional
maxima on the other two curves), and the radial distribu-
tion function of the = 0n level begins to look qualitatively
like the radial distribution function of the = 1n − level with
two maxima.
Further, Fig.5(b) implies that the peak in the radial dis-
tribution function of the = 0n level near the impurity is
redistributed among the = 0m states of the = 1, 2, ...n − −
Landau levels. Obviously, this is an analog of the phenom-
enon of the diving into continuum for a supercritical
charge in the absence of a magnetic field. In the latter case,
the lowest bound state dives into the lower continuum pro-
ducing a resonance whose wave function can be consid-
ered as redistributed over the lower continuum states with
energies of the order of the resonance width γ . All wave
functions from this region have an additional sharp peak
near the origin. As we see, when magnetic field is present,
there is a similar redistribution of the profiles of radial dis-
tribution functions near the impurity (note that as the im-
purity charge increases, the “redistribution” region shifts
Fig. 5. (Color online) The radial functions of the electron density of the = 0m state for the Landau levels = 0n (1), = 1n − (2), and
= 2n − (3) and three different values of the impurity charge: = 0.7ζ (a), = 1.3ζ (b), = 1.9ζ (c).
502 Low Temperature Physics/Fizika Nizkikh Temperatur, 2018, v. 44, No. 5
Electron states in the field of charged impurities in two-dimensional Dirac systems
down to the lower Landau levels). According to Fig. 5(c),
the blue curve (1) representing the electron density is now
similar to the red dashed curve (2) in Fig. 5(a) and the red
dashed curve (2) is similar to the green dot-dashed curve (3)
in Fig. 5(a).
So far we did not take into account the screening of
a charged impurity due the polarization effects in graphene
to which we turn our attention in the next subsection.
3.2. Tuning the screening of charged impurity
with chemical potential
Experimentally, as shown in Sec. 3.3, the strength of
a charged impurity and splitting of Landau sublevels with
different orbital momenta in a magnetic field can be very
effectively tuned by a gate voltage [36]. In this subsection,
we theoretically study this phenomenon by taking into ac-
count the polarization in a magnetic field which is con-
trolled by the chemical potential due to gate voltage.
It is natural to attribute the variation in the strength of
the impurity potential to the screening properties of the 2D
electron system. To describe this effect theoretically, we
follow Ref. 67 and use the polarization function calculated
in the absence of a charged impurity. Then the correspond-
ing Poisson equation, which defines the screened impurity
potential, reads
2
(2)
2 tot
2( ) = ( )D
ZeV x xπ
−∆ − δ −
κ
2
2
tot
2 ( ; ) ( ),e d Vπ
− Π − µ
κ ∫ y x y y (56)
where ( ; )Π − µx y is the static polarization function calcu-
lated by using the wave functions of free electrons in a mag-
netic field. Notice the presence of the pseudodifferential
operator 2D−∆ in the equation above, which is necessary
to correctly describe the Coulomb interaction in a dimen-
sionally reduced electrodynamic system [70–73].
Since Eq. (56) is algebraic in momentum space
2 2
tot
2 2| | ( ; ) ( ) = ,e ZeV
π π
+ Π µ − κ κ
q q q (57)
we easily find the screened impurity potential in coordinate
space
2 2
tot 2
exp( )( ) = =
2 2| | ( ; )
Ze d q iV
e
−
κ π π
+ Π µ
κ
∫
qrx
q q
2
0
2
0
( | |)
= ,
2 ( ; )
q J qZe dq
eq q
+∞
−
κ π
+ Π µ
κ
∫
x (58)
where = | |q q . The static polarization function at zero
temperature has the form [74],
2 2
2
=0 =
( ; ) = ( )
24
ncf B
nn n
nB
N q lq Q M
l
λλ
Γ
λ ±
Π µ δ µ −λ − π
∑∑
2 2
, =0 , =
( ) ( )
,
2
nc
n nB
nn
n nn n
n n
M Mq lQ
M M
′ ′λλ Γ Γ
′
′′ ′λ λ ±
′ ′λ ≠λ
′θ µ −λ − θ µ −λ − ′λ − λ
∑ ∑
(59)
where = 2 /n F BM n lv are the Landau level energies,
and we introduced the ultraviolet cutoff cn because of
the divergence of the sum over the Landau levels. Since
the bandwidth is finite in graphene, cn is estimated as
4= 10 / [ ]cn B T [75,76]. As in experiment [36], we consider
the system of two superposed graphene layers twisted away
from Bernal stacking by a large angle. This does not affect
the spectrum of single-layer graphene but results in an ad-
ditional twofold layer degeneracy: the factor = 2 2 = 4f s lN
takes into account spin degeneracy and the presence of the
second graphene layer. In experiment, this setup makes
possible to diminish the random potential fluctuations due
to substrate imperfections. The smeared delta function
2 2( ) = ( / ) / ( )x xΓδ Γ π + Γ and the step function
1 1( ) = ar ctg
2
xxΓ
θ + π Γ
account for the finite width of Landau levels, and the func-
tions ( )nnQ y′λλ
′ are defined as
___________________________________________________
2
0, < | | | || | > <
0, 1<< <> >
(1 ) ! ( 1)!
( ) = e ( ) (1 ) ( ) ,
! ( 1)!
n n n n ny n n
nn nn n
n nQ y y L y L y
n n
′ ′′ ′ − −λλ − −
′ −
′+ λλ δ − ′+ λλ − δ
−
(60)
______________________________________________
where < = min( , )n n n′ , > = max( , )n n n′ and ( )m
nL y are the
generalized Laguerre polynomials. The first term in Eq. (59)
describes the contribution from the intralevel transitions
while the second term represents the contribution from the
interlevel transitions. For a small width of Landau levels
the first term looks like a sequence of delta functions and
contributes only when the chemical potential lies inside
Landau levels. At small wave vectors ( 1
Bq l−
) the polari-
zation function (59) behaves as [74]
2
2( ; ) ( ),
2 e
TFq q dqκ
Π µ +
π
(61)
Low Temperature Physics/Fizika Nizkikh Temperatur, 2018, v. 44, No. 5 503
E.V. Gorbar, V.P. Gusynin, and O.O. Sobol
where
2
02
=0 =
= (2 ) ( )
ncf
TF n n
n
e N
q M
l
Γ
λ ±
− δ δ µ −λ
κ
∑∑ (62)
is the Thomas–Fermi wave vector which determines the
strength of the long-wavelength screening, and parameter
d is given by
2
0
=0 =
= (4 ) ( )
2
ncf
n n
n
e N
d n MΓ
λ ±
− + δ δ µ −λ −
κ ∑∑
2 1
1
3
=0 , =
( ) ( )
.
2 2 ( 1 )
ncf n n
F n
e N l M M
n n
−
Γ + Γ
′λ λ ±
′θ µ −λ − θ µ −λ
−
κ ′λ + −λ
∑ ∑
v
(63)
In fact, the static polarization function (0, )Π µ is pro-
portional to the density of states which at finite scattering
rate has the form of a series of broaden Landau levels [74].
Figure 6 illustrates the dependence of the static polarization
function (59) and its two leading long-wavelength terms
(62) and (63) on the chemical potential. We plot for compari-
son the unscreened potential and the screened potential (58)
of the impurity in Fig. 7. Let us consider the case where
the chemical potential is situated between Landau levels.
Then the Thomas–Fermi wave vector (62) is close to zero
(Figs. 6(a), (b)) and the Coulomb potential of the impurity
is weakly screened, although even in this case graphene
contributes to the total dielectric function at large and inter-
mediate momenta, which effectively diminishes the charge
of the impurity and the screened potential. Indeed, while
the screened potential tends to its bare value at r →∞ , it
is weakened for small and intermediate distances (see the
line (2) in Fig. 7). On the other hand, when the chemical
potential lies inside any given Landau level, the screening
is much more effective due to large TFq (see the line (3) in
Fig. 7) providing an excellent means of controlling the ef-
fective charge of impurity by the gate voltage which is di-
rectly related to the chemical potential µ. Moreover, the
coefficient d in Eq. (63) in this case is negative (see
Fig. 6(c)) which means that ( ; )qΠ µ has a nonmonotonic
momentum dependence with a peak at = 0q . This behavior
of the polarization function leads to the oscillations of the
screened potential (line (3) in Fig. 7) with the sign change
(i.e., the overscreening of the Coulomb potential) at inter-
mediate distances of the order of several magnetic lengths.
In Ref. 67 the backreaction of the charged impurity on
the polarization properties was also taken into account.
Although the qualitative picture of screening is the same, it
was shown that due to the downshift of the energy levels,
the polarization function no longer remains symmetric with
respect to the exchange µ ↔ −µ. These features of a charged
impurity in graphene in the magnetic field are clearly ob-
served in the recent experiments [36,37].
It should be noted that the approximation of noninter-
acting electrons may become invalid when the chemical
potential lies inside the Landau level. Indeed the electron-
electron interactions could lead in sufficiently clean gra-
phene specimen to such interesting phenomena as the frac-
tional quantum Hall effect. Then, the chemical potential
Fig. 6. (Color online) The dimensionless polarization function = (4 / ) ( ; )F B fl N qΠ π Π µ
v as a function of the chemical potential and the
wave vector (a) and the two coefficients 2= (2 / e )TF F B f TFq l N qκ
v and 2= (2 / e )F f Bd N l dκ
v of its expansion (61) at small wave
vectors (b), (c).
Fig. 7. The unscreened regularized Coulomb potential (1) and
the screened potential of the impurity as a function of / Br l in
the cases where the chemical potential is situated between Lan-
dau levels (2) and lies inside the zeroth Landau level (3).
504 Low Temperature Physics/Fizika Nizkikh Temperatur, 2018, v. 44, No. 5
Electron states in the field of charged impurities in two-dimensional Dirac systems
cannot be tuned continuously and instead jumps from one
plateau to another. Although our analysis becomes inappli-
cable in the fractional quantum Hall regime, the experi-
mental results in [36] show that the conclusion about the
maximal screening remains unchanged.
3.3. Screening charged impurities and lifting the orbital
degeneracy in graphene by populating Landau levels
Charged impurities in undoped gapless graphene pro-
duce a spatially localized signature in the density of states
(DOS) which is readily observed with scanning tunneling
microscope and spectroscopy (STM+STS)[77]. This effect
is especially important in the presence of a magnetic field
when the quantization of the 2D electronic spectrum into
highly degenerate Landau levels (LL) gives rise to the quan-
tum Hall effect. In this regime charged impurities are expect-
ed to lift the orbital degeneracy causing each LL in their
immediate vicinity to split into discrete sublevels [62,78].
By making use of high quality gated graphene devices
in a magnetic field, it was shown in Ref. 36 that the
strength of a charge impurity, as measured by its effect on
the electron spectrum, can be effectively controlled by tun-
ing the LL occupation with a back gate voltage. The LL
spectra were obtained by measuring the bias voltage depend-
ence of the differential tunneling conductance, /dI dV ,
which is proportional to the DOS, ( , )D E r , at the tip posi-
tion r . Here = ( ) /FV E E e− is the bias voltage and E is
the energy measured relative to the Fermi level, FE . For
almost empty LLs, the impurity is screened and essentially
invisible whereas at full LL occupancy screening is very
weak and the potential due to the impurity attains maxi-
mum strength. The underlying discrete quantum-mecha-
nical spectrum arising from lifting the orbital degeneracy
was experimentally resolved in the unscreened regime.
To explore the influence of the impurity on the LLs the
spatial evolution of spectra along a trajectory traversing it
for a series of gate voltages was studied. For certain gate
voltages the spectra become significantly distorted close to
the impurity, with the = 0n level (and to a lesser extent
higher order levels) shifting downwards toward negative
energies. The downshift indicates an attractive potential
produced by a positively charged impurity. Its strength, as
measured by the distortion of the = 0n LL, reveals a sur-
prisingly strong dependence on LL filling. In the range of
gate voltages 15 V < < 9 VgV− corresponding to filling
the = 0n LL the distortion grows monotonically with fill-
ing. At small filling the distortion is almost absent indicat-
ing that the impurity is effectively screened attaining its
maximum value close to full occupancy. At full occupancy
the = 0n level shifts by as much as 0.1 eV≈ indicating
that the effect would survive at room temperature. The
spectral distortion is only present in the immediate vicinity
of the impurity. Farther away no distortion is observed for
all studied carrier densities.
The variation of the impurity strength with filling is re-
lated to the screening properties of the electron system.
They were studied theoretically in Sec. 3.2. For a positively
(negatively) charged impurity and almost empty (full) LLs,
unoccupied states necessary for virtual electron transitions
are readily available in the vicinity of the impurity, result-
ing in substantial screening. By contrast for almost filled
(empty) LLs, unoccupied states are scarce, which renders
local screening inefficient.
The downshift is the largest for the = 0n and = 1/ 2j −
state and diminishes with increasing | |j and/or n. The
local tunneling DOS was calculated in Ref. 36 assuming
a finite linewidth Γ
†( , ) = 4 ( ) ( ) ( ),nj njnj
nj
D E E EΓδ − ψ ψ∑r r r (64)
where 2 2( ) = / [ (( ) )]nj njE E E EΓδ − Γ π − + Γ represents a
broadened LL. The peak intensity is determined by the pro-
bability density † ( ) ( )njnjψ ψr r and is position dependent. If
< njEΓ ∆ ( njE∆ defines spacing between adjacent levels),
the discreteness of the spectrum is resolved, but for njEΓ ≥ ∆
the peaks of adjacent states overlap and merge into a con-
tinuous band.
Thus, even if the spectrum is discrete, but the resolution
insufficient or if impurities are too close to each other, the
measured ( , )D E r will still display “bent” LLs, whose en-
ergies seemingly adjust to the local potential. In particular,
upon approaching the impurity the = 0n LL splits into
well resolved discrete peaks connected with specific orbital
states. As it could be seen from the experimental plots in
Ref. 36, the states 0 ( )jψ r with = 1/ 2, 3 / 2, 5 / 2j − − − are
well resolved close to the impurity but higher order states
are less affected and their contributions to ( , )D E r merge
into a continuous line. Similarly, the discreteness of the
spectrum is not resolved for 0n ≠ LLs that is consistent
with the weaker impurity effect at larger distances. For
partial filling ( = 5 V,0 VgV − ) as screening becomes more
efficient and orbital splitting is no longer observed, the
unresolved sublevels merge into continuous lines of “bent”
Landau levels. The ability to tune the strength of the impu-
rity in-situ demonstrated in Ref. 36 opens the door to ex-
ploring Coulomb criticality and to investigate a hitherto
inaccessible regime of criticality in the presence of a mag-
netic field [62,78].
The tuning of the effective charge due to the polariza-
tion effects was studied theoretically by the three of us in
Ref. 67. Numerically integrating the Dirac equation and
determining the energies of several first Landau levels in
the screened impurity potential (58), the local density of
states along the line cuts across the impurity was deter-
mined and is plotted in Fig. 8. The corresponding LDOS is
in good agreement with the experimental results of Ref. 36
(see Fig. 3 therein).
Low Temperature Physics/Fizika Nizkikh Temperatur, 2018, v. 44, No. 5 505
E.V. Gorbar, V.P. Gusynin, and O.O. Sobol
4. Two-electron bound states near a Coulomb impurity
In the weak interaction regime 0.4α ≤ , it was found in
Ref. 79 that a pair of repulsively interacting Dirac fermions
in graphene in the attractive potential of a Coulomb impu-
rity with charge Ze forms a two-body bound state localized
near the impurity. It could be observed by means of STM
experiments similar to those previously reporting super-
critical behavior in graphene [16,36,51] and trapped elec-
tron states in electrostatically defined graphene dots [80,81].
The negatively charged two-electron hydrogen ion H −
represents a classic problem of nonrelativistic quantum
mechanics [82–87]. As it was shown in Ref. 84, there ex-
ists a single bound state in three spatial dimensions. Chan-
drasekhar proposed to construct a trial wave function for
the ground state of H − as follows [82]:
1 2 1 21 2( , ) = e e ,ar br br arr r − − − −Ψ − ε (65)
where = | |l lr r , = 1,2l denote the distance of the corre-
sponding electron to the nucleus, a and b are the variational
parameters, and = 1ε corresponds to a spin singlet/triplet
state, respectively. The variational calculation shows that
the minimal energy for a two-body bound state is obtained
for a b≠ in the spin singlet configuration ( = 1ε − ).
The nonrelativistic 2D counterpart of the above system,
which is the D− problem, describes a donor impurity ion
with two electrons in a 2D semiconductor quantum well
[88–93]. The effects of quantum confinement on two-body
bound-state energies have been studied experimentally in
Ref. 94. In the absence of a magnetic field, there exists
only a single bound state in the spin singlet sector. The D−
problem is also similar to the negatively charged exciton
(X − ) problem, which was experimentally studied in quan-
tum wells [95].
The corresponding 2D relativistic problem could be re-
alized in gapped graphene monolayers (or topological insu-
lator surfaces) with a Coulomb impurity. However, it was
found long ago that the relativistic H − problem is subtle
because the single-particle Dirac Hamiltonian is unbound-
ed from below [96–98]. In order to set a physically and
mathematically well-posed problem, it is necessary to pro-
ject the interaction Hamiltonian onto the states with positive
energy. It can be devised for interacting Dirac fermions in
graphene if (i) a single-particle gap exists ( > 0∆ ), and (ii)
electron-electron interactions are weak, see Refs. 99–102.
We follow the derivation in Ref. 79 and consider the in-
teracting two-particle problem for a gapped graphene mo-
nolayer in the presence of a charged impurity. The cor-
responding Dirac–Coulomb Hamiltonian at a given valley
reads
2
=1,2
= ( ) .b
l
H H l V+∑ (66)
Here ( )H l is the single-particle Dirac Hamiltonian (1) for
particle = 1,2l with the single-particle potential of charge
Z impurity at the origin
2
= = g FZZeV
r r
α
− −
κ κ
v
(67)
and 2bV is the standard two-body Coulomb interaction
which equals
2
2
1 2
= .
| |b
eV
κ −r r
(68)
However, as was mentioned above, the bound state problem
is not well posed for the Dirac–Coulomb Hamiltonian of
two particles (66) because the spectrum of the single-
particle Dirac Hamiltonian is unbounded from below. A si-
milar problem occurs in the study of relativistic effects in
Fig. 8. (Color online) The local density of states calculated numerically in Ref. 67 is plotted at four values of gate voltage along the line
cuts across the impurity.
506 Low Temperature Physics/Fizika Nizkikh Temperatur, 2018, v. 44, No. 5
Electron states in the field of charged impurities in two-dimensional Dirac systems
the helium atom as was found long ago by Brown and Ra-
venhall [96]. According to Sucher [99,100], the Hamilto-
nian H in Eq. (66) has to be replaced by the projected
Hamiltonian [101,102]
2= (1) (2) ,bH H H V+ + ++ + Λ Λ (69)
with the projection operator = (1) (2)+ + +Λ Λ Λ . Since
2 1/2( ) = [ ( )]l H l is positive definite, the single-particle
operator [ ]( ) = ( ) ( ) / 2 ( )l l H l l+Λ + , obviously, projects
onto the positive energy states. As shown in Refs. 99–102,
the projected Hamiltonian H+ takes into account the most
important effects of the electron-electron interaction. Fur-
ther, due to the presence of a band gap, the replacement
H H+→ is reliable for the ground state of the system in
the limit of weak Coulomb repulsion. Moreover, the pro-
jection guarantees that the Hamiltonian H+ can possess the
genuine two-particle bound states.
Taking the two-particle wave function in a factorized
form tot = |Φ Φ χ〉, where | χ〉 is the normalized spin part
(singlet or triplet), its spatial part equals [79]
1 2 1 2 1 2( , ) = ( ) ( ) ( ) ( ),I O O IΦ Ψ Ψ −εΨ Ψr r r r r r (70)
where = 1ε for the spin singlet/triplet sector and IΨ and
OΨ are the normalized ground-state eigenspinors of the
2D relativistic hydrogen problem with charges Z replaced
by variational parameters IZ and OZ . These eigenspinors
are given by
( )( 1)/2 1
( ) = 2 e .
(1 ) e 1
p r
i
p
r p r
i
λγ − −λ λ λ
λ λ φ
λ λ
+ γ
Ψ
πΓ + γ − γ
(71)
Here = ,I Oλ and
2 2= 1 4 ,Zλ λγ − α = 2 / ( ).Fp Zλ λ∆ α v (72)
For the energy functional
tot tot
tot tot
| | | |
( , ) = = ,
| |I O
H HE Z Z + +
ε
〈Φ Φ 〉 〈Φ Φ〉
〈Φ Φ 〉 〈Φ Φ〉
(73)
calculating the matrix elements explicitly, we obtain the fol-
lowing energy functional (for details see Ref. 79):
dir exc
2 2
2 2
( )( )
( , ) = ,
1 1
b b
I O
Z Z V SU V V
E Z Z
S S
λ λ
ε λ
λ
− − ε − ε
∆γ + +
− ε − ε
∑
(74)
where the overlap integral equals
(1 )(1 ) (1 )(1 )
=
2
I O I OS
+ γ + γ + − γ − γ
×
(1 )/2(1 )/2
1 ( )/2
(1 ( ) / 2) (2 ) (2 )
,
(1 ) (1 ) ( )
OII O I O
I OI O I O
p p
p p
+γ+γ
+ γ +γ
Γ + γ + γ
×
Γ + γ Γ + γ +
(75)
one-particle matrix elements are
1= | ( / ) | = 2 ,
p
V r λ
λ λ λ
λ
〈Ψ α Ψ 〉 α
γ
(76)
1= | ( / ) | = 2 ,I O
I O
I O
p p
U r S
+
〈Ψ α Ψ 〉 α
γ + γ
(77)
and two-particle matrix elements
22dir
2 1 2 1 2
12
= ( ) ( ) ,b I OV dr dr r r
r
α
Ψ Ψ∫ (78)
† †exc
2 1 2 1 1 2 2
12
= ( ) ( ) ( ) ( )b O II OV dr dr r r r r
r
α Ψ Ψ Ψ Ψ ∫ (79)
could be represented in terms of elliptic functions (see
Ref. 79).
It should be noted that the energy functional (74) in-
cludes the matrix elements of the full interaction operator
rather than those of the projected operator, 12( / )r+ +Λ α Λ ,
which are more difficult to obtain and would require a de-
tailed numerical analysis. Both matrix elements coincide if
the trial wave function has vanishing projection onto the
negative energy eigenfunctions of DH . In fact, in Ref. 79,
it was verified that for 0.4α the cumulative weight of
negative energy states in the trial wave function is very
small ( 1% ). Indeed, negative energy states will only be
important if typical interaction matrix elements can over-
come the band gap 2∆. For small α, one therefore expects
at most small quantitative corrections in the bound-state
energy due to this approximation.
The energy functional (74) possesses the following fea-
tures. First of all, it is symmetric under an exchange of its
arguments. Second, for the spin-singlet case ( = 1ε − ), this
energy is bounded from below for < 1/ 2Zα . Third, for
small α, ( , )I OE Z Zε reduces to the corresponding nonrela-
tivistic energy functional for the D− problem in 2D semi-
conductors [92]. However, in contrast to the nonrelativistic
case, ( , )I OE Z Zε is not homogeneous in α, and hence the
bound-state energy explicitly depends on α. As in the non-
relativistic case, this energy minimum is realized for une-
qual values of IZ and OZ .
It was shown in Ref. 79 that the energy functional for
the singlet state with = 1Z has a minimum located below
the threshold, i.e., the binding energy is positive. In addi-
tion, there exists a two-body bound state. The situation is
different in the spin triplet sector, where the variational
approach predicts that the energy functional has a mini-
mum whose energy is above the threshold and, thus, does
not describe a bound state. For = 1Z , the minimum is at
< 1OZ and > 1IZ (or vice versa, due to the symmetry of
Eε). Physically, one quasiparticle partially screens the im-
purity charge seen by the other quasiparticle.
Also in Ref. 79 the authors calculated the probability
density and the pair distribution function for the bound
state, focusing on the two-body spin singlet state. They sug-
gest that the bound state can be accessed experimentally,
e.g., by means of STM techniques.
Low Temperature Physics/Fizika Nizkikh Temperatur, 2018, v. 44, No. 5 507
E.V. Gorbar, V.P. Gusynin, and O.O. Sobol
5. Two-center problem
Although, naively, the supercritical instability should be
easily realized for charged impurities with > 1Z , its exper-
imental observation remained elusive due to the difficulty
of producing highly charged impurities. However, one can
reach the supercritical regime by collecting a large enough
number of charged impurities in a certain region. As we
saw in Sec. 2.3, this approach was successfully realized by
creating artificial nuclei (clusters of charged calcium di-
mers) on graphene [16] using the tip of a scanning tunnel-
ing microscope. It is ironic that in spite of much larger
value of coupling constant in graphene than in QED the
first observation of the supercritical instability in graphene
still required the creation of supercritical potentials from
subcritical charges like in the case of heavy nuclei colli-
sions in QED discussed in the Introduction. What crucially
differs the graphene experiments [16] from those in QED
is that the supercritical electric fields created by placing
together ionized Ca impurities are static unlike the fields
created in heavy nuclei collisions in QED. This makes
possible to observe and analyse reliably the supercritical
regime.
The Hamiltonian of the two-center problem is the same
as Hamiltonian (1) with the potential
2
1 2
1 2
( ) = ,
Z ZeV
r r
− + κ
r (80)
where 1,2 =| / 2 |r ±r R are the distances from electron to
impurities with charges 1,2Z . Since the experiments in
Ref. 16 were performed for impurities of the same type, we
will study in what follows the symmetric problem, i.e.,
1 2= =Z Z Z . The alignment of the charges with respect to
the origin = 0r is arbitrary due to translational and rota-
tional invariance of the free gapped Dirac Hamiltonian (1),
and we choose them located at ( / 2,0)R with R being
a distance between two charges.
The main difficulty in solving the Dirac equation with
two Coulomb centers in QED is that variables in this prob-
lem are not separable in any known orthogonal coordinate
system [103]. Unfortunately, this is true also for the Dirac
equation for two Coulomb centers in the (2 1)+ -dimen-
sional problem in graphene. Therefore, we apply the ap-
proximate methods such as the linear combination of atom-
ic orbitals (LCAO) technique and variational method.
5.1. LCAO approach for symmetric two-center problem
The LCAO method is well known and widely used in
molecular physics [104]. Wave functions in this method
are chosen as linear combinations of basis functions, where
the latter are usually the electron functions centered on the
corresponding atoms of the molecule. By minimizing the
total energy of the system, the coefficients of the linear
combinations are then determined. The LCAO approach
for the symmetric two-center Dirac problem in 2D was
applied in Ref. 105 (see also Ref. 106), where only the lowest
single-impurity bound state near each center is retained.
This approximation is expected to yield accurate ground-
state energies for large R [107,108], where the molecular
ground state is well approximated in terms of atomic orbitals.
In addition, as we show below, the exact result for 0R →
is also captured by the LCAO solution.
For a single impurity of charge Z , the lowest bound state
has the energy γ∆ with 2 2= 1 4Zγ − α . In the absence of
short-distance regularization, the supercritical threshold is
reached at = 1/ 2cZ α [21], therefore, < cZ Z is assumed
henceforth. The corresponding normalized spinor, which is
an eigenstate of the total angular momentum operator
= / 2z zJ i φ− ∂ + σ with eigenvalue 1/ 2, has the same
form as (71) and reads
( 1)/2
2 /
0
12 4( , ) = e ,
(1 ) e 1
Z r R
i
Z Z rr
RR i
γ−
− α ∆
φ
∆∆
+ γ α α Ψ φ πΓ + γ − γ
(81)
where = /FR∆ ∆v .
It is convenient to rewrite the interaction potential as
follows:
2
0 eff
1 2
1 1= ( ) ,eH H Z Z
r r
− + δ + κ
(82)
where eff=Z Z Zδ − and the effective charge effZ is intro-
duced. It is a variational parameter, which could be deter-
mined by minimizing the total energy.
Following the standard LCAO approach [104], we seek
the electron wave function |Φ〉 in terms of atomic orbitals,
|1〉 and | 2〉 centered near the Coulomb impurity at
( / 2,0)R , which depend on 1r and 2r , respectively, i.e.,
1 2| = |1 | 2v vΦ〉 〉 + 〉 . The atomic orbitals are chosen as
eigenstates (81) in the field of a single impurity of charge
effZ . The Dirac equation is thereby reduced to a linear sys-
tem of equations for 1v and 2v , and the energy eff= ( )E E Z
follows from the condition
11 12
21 22
det = 0,
H E H SE
H SE H E
− −
− −
(83)
where | |ijH i H j= 〈 〉 and the overlap integral
= 1| 2 = 2 |1S 〈 〉 〈 〉 .
Defining the Coulomb integral,
1 1
2 1= 1| |1 = 2 | | 2 ,C r r− −〈 〉 〈 〉 (84)
and the resonance integral,
1 1
1,2 1,2= 1| | 2 = 2 | |1 ,A r r− −〈 〉 〈 〉 (85)
all matrix elements in Eq. (83) can be written in compact
form:
508 Low Temperature Physics/Fizika Nizkikh Temperatur, 2018, v. 44, No. 5
Electron states in the field of charged impurities in two-dimensional Dirac systems
2
2
11 22 eff
2
12 21
= = 4 / ,
e= = ( ) .
eH H Z Z Z C
H H S Z Z A
γ∆ − δ α ∆ γ −
κ
γ ∆ − + δ
κ
(86)
While S , C , and A can be directly evaluated [107,108] in
the 3D Dirac problem, the 2D case is, unfortunately, more
involved. In order to compute them, it is useful to employ
elliptic coordinates [46]
1 2= [1, ),
r r
R
+
ξ ∈ ∞ 1 2= [ 1,1].
r r
R
−
η ∈ − (87)
In these coordinates the integrals could be written as fol-
lows:
1
=
(1 )
uS
γ+
×
πΓ + γ
1 2 2 ( 1)/2
2 2
2 2
1 1
( ) e 1 (1 ) ,
( 1)(1 )
u
d d
∞ γ− − ξ
−
ξ −η × ξ η ξ − + γ −η
ξ − −η
∫ ∫ (88)
12=
(1 )
uA
R
γ+
×
π Γ + γ
1 2 2 ( 3)/2
2 2
2 2
1 1
( ) e 1 (1 ) ,
( 1)(1 )
u
d d
∞ γ− − ξ
−
ξ ξ −η × ξ η ξ − + γ −η
ξ − −η
∫ ∫ (89)
11 ( )
2 2
1 1
2 ( ) e= ,
(1 ) ( 1)(1 )
uuC d d
R
∞γ+ γ − ξ+η
−
ξ + η
ξ η
π Γ + γ ξ − −η
∫ ∫ (90)
where eff= 2 /u RZ R∆α and 2
eff= 1 (2 )Zγ − α .
Using the above integral representations for S, C, and A
(or the corresponding series from Ref. 105), it is numeri-
cally straightforward to obtain an LCAO estimate for the
ground-state energy eff( )E Z for given effZ . Then the min-
imal energy is determined, realized for *
eff =Z Z , where
the numerical search is aided by noting that ( )E Z depends
quadratically on *
effZ Z− . Numerical results in Ref. 105
show that the LCAO ground-state energy, ( )E R , matches
the expected single-impurity values γ∆ in both limits, na-
mely (i) for R →∞ with impurity charge Z, where we have
two decoupled copies of the single-impurity problem, and
(ii) for 0R → , where both centers conspire to form a single
Coulomb impurity of charge 2Z . Furthermore, it was shown
that the optimal effective charge *Z nicely matches both
limits as well.
Choosing larger Z such that = / = 2cZ Z Zζ α is within
the bounds 1/ 2 < < 1ζ , the supercritical regime can be real-
ized by decreasing R through a transition value, cr=R R .
At the critical distance, the ground-state energy reaches the
Dirac sea, cr( ) = ,E R −∆ and for cr<R R , the two-center
system with subcritical individual impurity charge be-
comes supercritical. The LCAO results in Ref. 105 show
that, in practice, Z has to be chosen quite close to = 1 2cZ / α,
since otherwise crR becomes extremely small. This conclu-
sion seems also in agreement with the reported experi-
mental observations of supercriticality [16,36], where dif-
ferent ions first had to be pushed closely together, thereby
forming charged clusters, before supercriticality appears.
5.2. Variational method
Another means to study the supercritical instability of
two Coulomb centers in graphene is the variational method
which was applied to the corresponding two centers prob-
lem in QED in Ref. 103. As noted in Ref. 109, in order to
obtain a satisfactory accuracy it is necessary that trial func-
tions correctly reproduce the asymptotics of the exact solu-
tion at infinity and near the charged impurities. These
asymptotics could be found from the direct analysis of the
Dirac equation ˆ =H EΨ Ψ . For two-component spinor
( ) = ( , )TΨ φ χr , expressing χ in terms of φ leads to the
following second order equation for the φ component of
the Dirac spinor:
2 2
2 2
2
( )( ) = 0.
( )
x y
F
V Vi
E Vx y i
E V x y
∂ ∂
−
∂φ ∂φ − − ∆∂ ∂
∂ + ∂ φ+ + + φ − + ∆ ∂ ∂ v
(91)
According to Refs. 3, 4, the supercritical instability takes
place when the bound state with the lowest energy dives
into the lower continuum. This occurs when =E −∆ . For
this solution, let us consider the asymptotic at large r R ,
where the potential equals
( )
2
23 5
1 1= (cos ) ,
4
F
RV P O
r r r
−ζ + ϕ +
r v (92)
= 2Zζ α is a dimensionless charge, and 2 ( )P x is the Legendre
polynomial ( )nP x with = 2n . In what follows we consider
the case when charges of impurities are subcritical whereas
their total charge exceeds a critical one, 1/ 2 < < 1ζ . The case
< 1/ 2ζ corresponds to the situation when the total charge is
less than a critical one and is not relevant for the supercri-
tical regime.
Neglecting the quadrupole and higher order multipole
terms in the potential (this corresponds to the monopole
approximation) Eq. (91) reduces to the following equation
for ( )rφ :
2
2
2 2 = 0,m
r rr
ζ ζ′′ ′φ + φ + − φ
(93)
where = / ( )Fm ∆ v . The decreasing at infinity solution is
expressed in terms of a Macdonald function,
1/2 2
1( ) = ( 8 ), = 4 1.ir C r K m r−
βφ ζ β ζ − (94)
Its asymptotic
3/4
asym 1( ) = exp ( 8 ),r C r m r r−φ − ζ → ∞ (95)
agrees, of course, with the asymptotical behavior of a solu-
tion to the Dirac equation of one center with charge 2Ze.
Low Temperature Physics/Fizika Nizkikh Temperatur, 2018, v. 44, No. 5 509
E.V. Gorbar, V.P. Gusynin, and O.O. Sobol
This shows also that the level which reached the boundary
of the lower continuum remains localized.
In order to find the asymptotic of the solution in the vi-
cinity of two Coulomb centers, it is conventional and con-
venient to use the elliptic coordinate system (87). We note
that for charges of impurities such as < 1/ 2Zα ( < 1ζ )
there is no “collapse” in the Coulomb field of one impurity
[12–14,21,110–112], therefore, it is not necessary to cut
off the potential at small r and the impurities may be con-
sidered as point-like. Therefore, for simplicity, in what fol-
lows, we will consider the nonregularized Coulomb poten-
tial. In elliptic coordinates, it has the form
( ) 2 2
2
= .
( )
FV
R
ζ ξ
−
ξ −η
r v (96)
To find the asymptotic of φ in the vicinity of impurities,
i.e., for small 2 2ξ −η , we seek φ in the form ( , ) = ( )φ ξ η φ µ ,
where 2 2 2
1 2= = 4 /r r Rµ ξ −η . Near the impurities 1 0r →
or 2 0r → , i.e., 1ξ → and 1η→ ± and, consequently, 0µ → ,
we obtain the following equation:
2 2
2 2
2 = 0,
4
d d
dd
φ φ ζ
+ + φ
µ µµ µ
(97)
whose regular solution at 0µ → is
/2
imp 2( ) = ,C −σφ µ µ 2= 1 1 .σ − − ζ (98)
This asymptotic describes the behavior of the wave
function at the impurities positions. Since at large distances
r R , the variable µ equals 2 24 /r Rµ , solution (94)
can be rewritten as follows:
1/4 1/4
1( ) = (2 ).iC K m R−
βφ µ µ ζ µ (99)
Matching solutions (98) and (99) at the point = 1µ we
can find an approximate estimate of the critical distance
cr ( )R ζ as a function of ζ . We obtain the following trans-
cendental equation:
2 (2 )
2 1 1 = 2 .
(2 )
i
i
K m R
m R
K m R
β
β
′ ζ
− ζ − ζ
ζ
(100)
For 1mR , i.e., when the distance between the impuri-
ties is much less than the Compton wavelength of
quasiparticles, Eq. (100) can be simplified using the as-
ymptotic of ( )iK zβ for 0z → . Then we obtain the follow-
ing analytical solution:
2
1
cr
1 2 11 2= exp arg (1 ) ,cotmR i−
− − ζ − − Γ + β
ζ β β
(101)
where ( )zΓ is the Euler gamma function. It is amazing that
Eq. (101) coincides with the corresponding solution found
in QED for scalar particles [113]. Equation (101) for
24 1 1ζ − can be written in more simple form
cr 2
1 2= exp .
4 1
mR
π −
ζ ζ −
(102)
We find that the deviation of crR given by Eq. (102)
from that determined by Eq. (100) is rather small up to
= 0.8ζ . A numerical calculation of crR given by these
equations is presented in Fig. 9 in comparison with crR de-
termined in more refined calculations using a variational
method.
Clearly, the approximation we used is rather crude be-
cause it matches only the asymptotics and, in particular, it
does not take into account at all the nonsphericity of the
potential of two impurities described by 2 (cos )P θ and
higher harmonics in potential (92).
To set up the variational problem, we note that the dif-
ferential equation (91) can be obtained as an extremum of
the following functional:
2
1 2
2
( )[ ] = ( ) | | ,
( )F
E VS E V i dxdy
x y
−
∂φ ∂φ − − ∆ φ − + ∆ + − φ
∂ ∂
∫
v
(103)
under the condition that the norm *=N dxdyΨ Ψ∫ is con-
served (this condition is important for obtaining the correct
boundary conditions). Introducing a new field 1/2= W −ψ φ,
where =W E V− + ∆ , the functional [ ]S φ can be repre-
sented in the form specific for nonrelativistic quantum me-
chanics
2 *[ ] = | |
2
VS i
W
∇ ψ ∇ψ + ×∇ψ ψ − ∫
2* 2( ) | | ,
2
Vi U dxdy
W
∇ − ψ ×∇ψ + − ε ψ
(104)
where = ij i ja b× εa b , 2 2 2 2= ( ) / (2 )FEε − ∆ v is the effec-
tive energy, and the effective potential U is given by
Fig. 9. (Color online) The dependence cr ( )mR ζ given by
Eqs. (100) (1) and (101) (2), and calculated by variational method
with = 1N (3) and = 2N (4).
510 Low Temperature Physics/Fizika Nizkikh Temperatur, 2018, v. 44, No. 5
Electron states in the field of charged impurities in two-dimensional Dirac systems
2 2
2 2
2 3 ( )= .
4 82( )F
EV V V VU
W W
− ∆ ∇
+ +
v
(105)
The second and third terms in functional (104) describe the
pseudospin-orbit coupling with the field = / 2V W−∇F ,
they do not contribute for the ground state wave function
which is real. Functional (104) is bounded from below, so
one is in position to apply to it the variational principle. In
what follows we are interested in the case where the bound
state with the lowest energy crosses the boundary of the
lower continuum, so we put =E −∆ i.e., = 0ε . Then =W V−
and the functional [ ]S ψ is simplified.
In QED, the Ritz and Kantorovich methods were em-
ployed in order to solve the variational problem and find
a critical distance crR (see a discussion in Sec. 3 in Ref. 114).
In the Ritz method, the sought function ψ is expanded over
a fixed set of basis functions ( , ) = ( , )n nnx y c x yψ ψ∑ ,
where nc are variational constants. In the Kantorovich
method, = ( ) ( )n nnc x yψ ψ∑ , where ( )n yψ are fixed func-
tions, while ( )nc x are variable functions. Obviously, the
variational problem reduces to a system of linear algebraic
equations for nc in the Ritz method and to a system of line-
ar ordinary differential equations for ( )nc x in the Kanto-
rovich method.
According to Eq. (98), near impurities φ depends only
on 2 2 2
1 2= = 4 /r r Rµ ξ −η . At the large distances, ,r →∞
the variable µ →∞ and the asymptotic of φ is given by
Eq. (99). Therefore, both asymptotics of ψ depend only on
µ. In order that variational ansatz for ψ give appropriate
results, it is essential to take into account correctly the be-
havior of the exact solution near the Coulomb centers and
at infinity. We choose the two variables ,µ ν so that the
function ( , )x yψ has a singularity only in µ. Then using
the following ansatz in the Kantorovich method:
1
=1
= ( ) ,
N
k
k
k
−ψ ψ µ ν∑ (106)
where ( )kψ µ are variable functions of µ and ( , )ν ξ η is
a fixed function of ξ and η, we can maximally correctly
take into account the behavior of the exact solution near
the Coulomb centers. Since a priori we do not know what
set of functions ( , )n x yψ is the best in the Ritz method, we
use like in the QED studies [103] the Kantorovich method.
Two choices of function ν were considered in QED [103]:
i) 2 2 2 2
1 2 1 2= / ( ) = ( ) / 4r r r rν η ξ −η − and ii) 2= =ν η
2 2
1 2= ( ) /r r R− . The obtained results were close. Here, we
will consider the case (i). Since the charges of impurities
are identical, 1 2= =Z Z Z , the wave function of the ground
state is symmetric under the inversion ,x x→ − y y→ − ,
therefore, the change of the variables ,x y to ,µ ν is per-
formed by means of the formulas
1/2= ( 1), = [( ( 1) 1)(1 )] .
2 2
R Rx yµ ν ν + µ ν + − −µν (107)
Inserting ansatz (106) in Eq. (104) and integrating over ν,
we obtain
( * *
, =1 0
( ) = 4
N
N kl k l kl k l
k l
S d P Q
∞
′′ψ µ ψ ψ + ψ ψ +∑ ∫
)†* * ,kl k l k lklR R ′′+ ψ ψ + ψ ψ (108)
where ,P Q, and R are N N× matrices which depend on µ
2 2
0
( ) = ( ) | | ( , ) ,k l
klP J f d
∞
+ −µ µ ν µ ν ν∫ ∇ (109)
42
0
( ) = ( ) ( 1)( 1) k l
klQ l k
∞
+ −
µ ν − − ν −
∫ ∇
3 2( ) 2 | | ( , ) ,
2
k l k lVi l k U J f d
V
+ − + − − − × ν ν + ν µ ν ν
∇
∇ (110)
3 2
0
( ) = ( 1)
2
k l k l
kl
VR l i
V
∞
+ − + − µ µ ν − ν + × µ ν × ∫
∇
∇ ∇ ∇
| | ( , ) .J f d× µ ν ν (111)
Here ( , ) = (1 )[ (1 ) ( ( 1) 1) ( 1)]f µ ν θ −µν θ −µ θ µ ν + − + θ µ −
and 2| | = / (16 ( 1)( 1)(1 ))J Rµ ν ν + µ +µν − −µν is a Jaco-
bian, ∇ is a gradient with respect to Cartesian coordinates,
U is the effective potential. The explicit expressions for
µ∇ , ν∇ , and U are given in Appendix A in Ref. 117.
Minima of functional (108) are given by solutions of
the following set of Euler–Lagrange equations:
† = 0.k k
kl k kl k klkl
d dd P R Q R
d d d
ψ ψ
+ ψ − ψ − µ µ µ
(112)
The boundary conditions for functions kψ follow from
the requirement that the norm of the function ψ be finite.
The differential equation (112) and these boundary condi-
tions define our boundary value problem.
In the simplest case = 1N , we have
= 0,d dP Q
d d
ψ
− ψ µ µ
(113)
where ( ) =P µ πµ and ( )Q µ is expressed through the com-
plete elliptic integrals of the first and second kind:
2( 1) ( 1) 1( ) = (1 )K( ) K
8 2
mRQ
π ζ − ζ θ µ −
µ − + θ −µ µ + + µ µ µ
23 (2 1)(1 )(1 ) E( ) K( )
8 (1 ) 8
ζ + +µ
+ θ −µ µ − µ +
µ +µ µ
2
2
3 1 ( 1)(1 ) 3 1( 1) E K .
8(1 ) 4
ζ − +µ + µ
+ θ µ − − + µ µ µµ
(114)
Low Temperature Physics/Fizika Nizkikh Temperatur, 2018, v. 44, No. 5 511
E.V. Gorbar, V.P. Gusynin, and O.O. Sobol
We seek a wave function of the ground state which could
be chosen real. Therefore, the function R11, which is com-
pletely imaginary, does not appear in Eq. (114). The dif-
ferential equation (113) determines the wave function of
the critical bound state that just dives into the lower con-
tinuum. Since the wave function of a bound state tends to
zero at infinity, this translates in our case to the condition
( ) 0ψ µ → as µ →∞. The asymptotic of the wave function
near the impurities (where 0µ → ) is given by Eq. (98).
This equation completes the set-up of our boundary value
problem which allows us to determine the critical distance
crR between the impurities as a function of ζ . Since the func-
tion ( )Q µ is given in terms of the complete elliptic inte-
grals of the first and second kind, the differential Eq. (113)
cannot be solved analytically. We solve this equation nume-
rically by using the shooting method and proceed as follows.
We fix the wave function and its first derivative at certain
small µ using Eq. (98). Then, we fix ζ and solve Eq. (113)
numerically for different mR (note that since the function
( )Q µ depends only on the product mR, parameters m and R
cannot be separately varied). The critical distance crR (for
a given m) is then determined as R such that the wave
function ( )ψ µ tends to zero at infinity. Repeating this pro-
cedure for different ζ , we find how the critical distance
between the impurities depends on ζ . The corresponding
dependence crmR on ζ is plotted in Fig. 9 (line 3).
The accuracy of computation can be improved taking
> 1N in sum (106). In this case one should solve a set of
second-order differential equations. Since the shooting me-
thod is not well suited for this purpose, it is better then to
follow the corresponding calculations in QED in Ref. 115
and reduce the set of Eqs. (112) to a matrix Riccati equa-
tion, which can be solved by the Runge–Kutta method. The
case = 2N was considered in Ref. 116 and the obtained
results are quite close to the case = 1N and are shown in
Fig. 9 (line 4).
5.3. Quasistationary states
When the distance between the impurities becomes
smaller than the critical one the bound state dives into the
lower continuum transforming into a resonance. The ener-
gy and width of this quasistationary state could be deter-
mined using the Wentzel–Kramers–Brillouin method in the
monopole approximation [117]. For distances >r R (or
more exactly r R ), the potential of the two-center prob-
lem is close to a spherically symmetrical one. Therefore,
we can consider one charged impurity with the charge 2Ze
and restrict our consideration only to the region r R≥ ,
where 1 is a dimensionless constant. This is known as
the monopole approximation.
For a spherically symmetric potential, the squared Dirac
equation could be rewritten in the Schrödinger-like form (33)
and the corresponding quasiclassical momentum ( , )k r E is
given by Eqs. (34)–(36). To find the energy of quasibound
states we use the Bohr–Sommerfeld quantization condition
0
/2 /2cr
( , ) = ( , ) = 2 ,
r r
R R
k r E dr k r dr n
− −
−∆ π∫ ∫
(115)
where = 1,2...n and r− is a boundary of the classically for-
bidden region determined by the equation ( , ) = 0k r E± and
0 = ( = )r r E− − −∆ .
For energies close to the boundary of the lower continu-
um, E → −∆ , we find
( , ) = ( , ),E R F Rζ −∆ ⋅ ζ
2 2 2 2
cr
2 2 2 2
1 2 1 2( , ) = 1 .
4 3 4 3
R j j j jF R
R
+ ζ − + ζ −
ζ + − + − ζ ζ ζ ζ
(116)
The width of quasistationary states apart from a preex-
ponential factor is determined by tunneling through the clas-
sically forbidden region. For energies close to the boundary
of the lower continuum, this gives the width
2
2 2
2 2exp 2 E j
E
Γ ∝ − π ζ − ζ −
− ∆
2 2
c
exp 2 ,
r
R j
R R
− π β − ζ − −
2 28 4 6 3= ,
24
j jζ + + +
β (117)
which tends to zero when E → −∆ or crR R→ .
6. The dipole problem and migration
of the wave function
The electron states for Dirac fermions in the field of
two oppositely charged ( Q± , =Q Ze) nuclei at distance R
(ZeR is the corresponding electric dipole moment) in
gapped graphene are described by Hamiltonian (1) with
potential
2
2 2 2 2
e 1 1( ) = ,
( / 2) ( / 2)
ZV
x R y x R y
−
κ + + − +
r (118)
where κ is the dielectric constant.
The corresponding Hamiltonian has an intrinsic parti-
cle-hole symmetry ˆ ˆ=H H+Ω Ω − , where the unitary opera-
tor = x xRΩ σ satisfies 2 = 1Ω ( xR is the operator of reflec-
tion x x→ − ). It follows then that an eigenstate ( , )E x yΨ
with energy E has a partner ( , ) = ( , ) =E Ex y x y−Ψ ΩΨ
= ( , )x E x yσ Ψ − with energy E− , hence, all solutions of
the Dirac equation come in pairs with E± .
This dipole problem was recently considered in
Refs. 105, 118 (the 3D Dirac equation with the electric di-
pole potential was also studied some time ago in Ref. 119).
It was shown that the point electric dipole potential (this
potential is defined as the limit 0R → of the finite-size
electric dipole potential (118) with fixed electric dipole mo-
512 Low Temperature Physics/Fizika Nizkikh Temperatur, 2018, v. 44, No. 5
Electron states in the field of charged impurities in two-dimensional Dirac systems
ment) accommodates towers of infinitely many bound states
exhibiting a universal Efimov-like scaling hierarchy and at
least one infinite tower of bound states exists for an arbi-
trary dipole strength. Notice that the Schrödinger equation
in two dimensions for the electron in the field of an electric
dipole also admits a bound state for any dipole strength [120]
unlike the three-dimensional case where a bound state exists
only when the dipole moment exceeds a certain critical
value (see, e.g., a discussion including a historical one in
Ref. 121). By combining analytical and numerical methods,
the authors of Ref. 118 found that the bound states do not
dive into the lower continuum because the positive and
negative energy levels first approach each other and then
go away. Actually this behavior is typical for an avoided
crossing [68], which forbids level crossing for two states
with the same quantum numbers. Since the bound states do
not dive into the lower continuum, the authors of Ref. 118
concluded that supercriticality is unlikely to occur in the
electric dipole problem in graphene.
We reconsidered the problem of supercriticality in our
paper [122] for the case of two oppositely charged impuri-
ties situated at finite distance (finite-size electric dipole).
By using the LCAO technique and variational Galerkin–
Kantorovich method, we showed that for sufficiently large
charges of impurities the wave function of the highest en-
ergy occupied bound state changes its localization from the
negatively charged impurity to the positive one as the dis-
tance between the impurities changes (both methods gave
similar results). The necessary condition for the instability
to occur is the crossing of the electron energy levels in the
field of single positively and negatively charged impurities.
This migration of the electron wave function of the super-
critical electric dipole is a generalization of the familiar phe-
nomenon of the atomic collapse of a single charged impurity
with holes emitted to infinity to the case where both elec-
trons and holes are spontaneously created from the vacuum
in bound states with two oppositely charged impurities thus
partially screening them.
6.1. Point electric dipole in graphene
Let us start from the case of 21/ r point electric dipole
potential, which was considered in Ref. 118. Far away from
the nuclei, r R , Eq. (118) is well approximated by the
point-like dipole form
2
2
cos( , ) = ( ) ,d F
pV r
r
θ
θ − v (119)
where 2 2 2= /( )Fp ZRe κ v is the quantity proportional to
the effective electric dipole moment. The 0r → singularity
requires regularization to avoid the usual fall-to-the-center
problem, see below. For nonrelativistic Schrödinger fermi-
ons, the dipole captures bound states only above a finite
critical dipole moment in three dimensions (3D) [119,120,
123–125]. However, a dipole binds states for arbitrarily
small p in the 2D Schrödinger case [120].
The corresponding 3D Dirac electric dipole problem
has not been discussed in (3+1)-dimensional QED presum-
ably because of the lack of heavy anti-nuclei preventing
its experimental realization in atomic physics. However, it
could be directly studied using STM spectroscopy in gra-
phene [16,36,51]. Bound states inside the gap, = ( )bE ± ∆ − ε
with binding energy bε ∆ , come in ( ,j ) towers of def-
inite “angular” quantum number, = 0, 1, 2,j , and parity
= ± (with 0j + ≥ ). The ( , )j tower is only present if
the dipole moment exceeds a critical value, ,> jp p , but
then contains infinitely many bound states. Since 0, = 0p + ,
there is at least one such tower. Bound states in the same
tower obey the scaling hierarchy
2 /, 1 ,
,
= e , = 1,2,
sb n j
b n
n
− π+ε
ε
(120)
where for p close to (but above) ,jp ,
,
,
2 , ( , ) = (0, ),
( )
( ) , > 0,j
j
p j
s p
p p j
∆ +
β − ∆
(121)
with 0.956β ≈ . Equation (120) agrees with the universal
Efimov law for the binding energies of three identical bos-
ons with short-ranged particle interactions [126–128]. Nu-
merical diagonalization of the Dirac equation in a finite
disc geometry indicates that as p increases, the bound
states approach = 0E without ever reaching it. The ab-
sence of zero modes was shown analytically in Ref. 118.
For energies close to the band edge = bE −∆ + ε with
| |bε ∆ , where > 0bε corresponds to bound states inside
the gap and for 2 2/ ( )Fp R ∆ v , the upper spinor com-
ponent stays always “small”,
1e
2
iF
ri
r
− θ
θ
φ ∂ + ∂ χ ∆
v
.
The equation, which determines the lower spinor com-
ponent, has the form
2
2( )
= 0,
2
F
bV
− ∇ − + ε χ ∆
v
(122)
with the 2D Laplacian 2∇ . We proceed with the potential
= dV V in Eq. (119), where Eq. (122) is solved by the ansatz
( , ) = ( ) ( )r F r Yχ θ θ . With separation constant γ , the angular
function satisfies an bε -independent Mathieu equation,
2
2 2 cos ( ) = 0,d p Y
d
+ γ − ∆ θ θ θ
(123)
which admits 2π-periodic solutions only for characteristic
values ,= ( )j pγ γ , where = ± is the parity, i.e.,
, ,( ) = ( )j jY Y−θ θ , and due to the anisotropy, = 0, 1, 2,j
differs from conventional angular momentum, with 0j + ≥ .
Using standard notation [46,129], the solutions to Eq. (123)
Low Temperature Physics/Fizika Nizkikh Temperatur, 2018, v. 44, No. 5 513
E.V. Gorbar, V.P. Gusynin, and O.O. Sobol
are expressed in terms of Mathieu functions 2ce j and 2se j ,
with eigenvalues 2 ja and 2 jb , respectively,
, 2 , 2
, 2 , 2
1( ) = ce ,4 , = (4 ),
2 4
1( ) = se ,4 , = (4 ).
2 4
j j j j
j j j j
Y p a p
Y p b p
+ +
− −
θ θ ∆ γ ∆
θ θ ∆ γ ∆
(124)
The characteristic values are ordered as
0, 1, 1, 2,< < < <+ − + −γ γ γ γ for given p. With ,= ( )j pγ γ ,
the radial equation reads
2
2 2 2
1 2 ( ) = 0.
( )
b
F
d d F r
r drdr r
∆εγ
+ − −
v
(125)
To regularize the fall-to-the-center singularity, the Dirichlet
condition 0( ) = 0F r is imposed at a short-distance scale
0r R≈ (on the level of the Dirac equation, this corresponds
to vanishing radial current at 0=r r ). It is found that this
regularization does not affect universal spectral properties
such as the Efimov law (120).
Let us now look for bound states, > 0bε . The solution
of Eq. (125) decaying for r →∞ is the Macdonald function
( )2 / ( )b FK rγ ∆ε v [46], and the condition 0( ) = 0F r
then yields an energy quantization condition within each
( ,j ) tower. Thereby the binding energies, , , , =b n jε
2 2 2
0= ( ) / (2 )n Fz r∆v , are expressed in terms of the posi-
tive zeroes, 1 2> > > 0z z , of
,
( )
j
K zγ
. Since only
( )isK z (with imaginary order) has zeroes [46], bound states
require , ( ) < 0j pγ . This condition is satisfied for ,> jp p
with
, ,( ) = 0.j jpγ (126)
The lowest few >0,jp resulting from Eq. (126) could be
found in Ref. 118. With increasing dipole moment, each
time that p hits a critical value ,jp , a new infinite tower
of bound states emerges from the continuum. Since
0, ( ) < 0p+γ for all p [129], we find 0, = 0p + : at least one
tower is always present. Explicit binding energies follow
from the small-z expansion of ( )isK z [46]. With the posi-
tive numbers , ,( ) = ( )j js p p−γ for ,> jp p , see
Eq. (121), we obtain
2 2 ( ) 2 /, ,
, , , 2
0
2
= e e ,
s n sj jF
b n j
r
ϕ − π
ε
∆
v (127)
where ( ) = (2 / ) arg (1 )s s isϕ Γ + . This becomes more and
more accurate as n increases. For n →∞ , in view of the
particle-hole symmetry, the energies accumulate near both
edges, , 0b nε → . Importantly, Eq. (127) implies the Efimov
scaling law announced in Eq. (120). This relation has its
origin in the large-distance behavior of the dipole potential,
and is thus expected to be independent of short-distance
regularization issues. A similar behavior has been predicted
for the quasi-stationary resonances of a supercritical Cou-
lomb impurity in graphene [12,21], and for 3D Schrödinger
fermions [123,125].
6.2. Migration of the wave function
in the finite dipole potential
In this subsection we consider the case of finite electric
dipole with potential (118). Since the variables in Dirac
equation with this potential are not separable in any ortho-
gonal coordinate system, we will utilize the LCAO meth-
od. It is convenient to work with dimensionless quantities
= /h H ∆ and = /E∆ε ∆ and use dimensionless coordi-
nates and distances defined in units of = /FR∆ ∆v . We
use also dimensionless coupling constant 2= e / ( )FZζ κv .
The LCAO method was used in Sec. 5.1 to solve the two
Coulomb centers problem in graphene. Here we apply LCAO
method to the dipole problem in graphene. As to the atomic
orbitals, we take the wave function of the lowest energy
bound state in the field of positively charged impurity and
the wave function of the highest energy bound state for
negatively charged impurity (these wave functions are re-
lated to each other by charge conjugation).
We begin our analysis with the Dirac equation for the
electron in graphene with one positively charged impurity
=p p ph ∆Ψ ε Ψ with the Hamiltonian
2 2
0
= ( ) .p x x y y zh i
r r
ζ
− σ ∂ + σ ∂ + σ −
+
(128)
(The Hamiltonian nh for the electron in the field of
negatively charged impurity is obtained from the Hamilto-
nian ph by the change of the sign of the last term in ph .)
We determine numerically the energy levels by using
the shooting method with regular boundary conditions at
= 0r for the wave functions and requiring that the wave
functions decrease at infinity.
The energy of the lowest (highest) electron bound state
with the total angular momentum = 1/ 2j in the regular-
ized Coulomb potential with the charge Ze+ ( Ze− ) is plot-
ted in Fig. 10 for different values of the regularization pa-
rameter 0r as a function of ζ . The levels which descend
from the upper continuum correspond to the positive charge
Ze+ while those which are pushed from the lower continu-
um and grow with ζ correspond to the negative charge
Ze− . These results are in accordance with calculations in
Ref. 21 (see Fig. 4 there) where slightly different regulari-
zation for the one Coulomb center potential was used which
admitted an analytical solution. They also reproduce quali-
tatively the behavior seen directly at the tight-binding level
on a honeycomb lattice [41]. For nonregularized Coulomb
potential with positive charge the lowest bound-state energy
is always positive, it reaches the value = 0ε for = 1/ 2ζ and
becomes purely imaginary for > 1/ 2ζ (the fall-to-center
514 Low Temperature Physics/Fizika Nizkikh Temperatur, 2018, v. 44, No. 5
Electron states in the field of charged impurities in two-dimensional Dirac systems
phenomenon [12–15]). For regularized Coulomb potential
with charge ( )Ze Ze+ − , the lowest (highest) bound-state
energy crosses = 0∆ε and dives into the lower (upper) con-
tinuum at certain value of the charge. For example, for
0 = 0.05r R∆, this happens for 1ζ ≈ .
Since the operator =c xU Kσ interchanges the ph and nh
Hamiltonians, =c p c nU h U h+ − , the electron levels in the
field of negatively charged center described by the Hamil-
tonian nh are obtained by the reflection ∆ ∆ε → −ε and
intersect with the levels of the Hamiltonian ph at = 0∆ε .
The corresponding critical value cζ when this happens will
play a crucial role in the behavior of energy levels in the
dipole potential because the behavior of these levels dra-
matically changes depending on whether < cζ ζ or > cζ ζ .
For chosen values of the regularization parameter 0r in
Fig. 10, the critical coupling 0= 0.6( = 0.01 )c r R∆ζ and
0= 0.7( = 0.05 )c r R∆ζ . In general, cζ increases with the in-
crease of 0r (see Fig. 1(b) and Eq. (31)).
We are ready now to consider the Dirac equation for
quasiparticles in graphene with two oppositely charged
impurities. The corresponding Hamiltonian has the form
2 2 2 2
0 0
= ( ) ,x x y y z
n p
h i
r r r r
ζ ζ
− σ ∂ + σ ∂ + σ + −
+ +
(129)
where 2 2
, = ( / 2)p nr x R y± + . We seek the wave function
as a linear combination (hybridization),
| = | | ,p p n nv vΨ〉 Ψ 〉 + Ψ 〉 (130)
of the wave functions pΨ and nΨ which are eigenstates of
the Hamiltonians ph and nh , respectively, with eigenvalues
0±ε , and 0ε is the energy of the lowest-energy electron bound
state in the field of one Coulomb center with the charge Ze+ .
Explicitly, the functions | ,|p nΨ 〉 Ψ 〉 with the total angular
momentum = 1/ 2j are given in polar coordinates by
( )
e ( )= , = ,
( )e ( )
ip n n
p ni p
np
f r g r
if ri g r
− θ
θ
Ψ Ψ −−
(131)
where , ,exp [ ] = ( / 2 ) /p n p ni x R iy r− θ ± − . The radial func-
tions ( )f r and ( )g r are computed numerically in the regu-
larized potential of one positively charged impurity.
It is crucial for our analysis below to use the electron
wave functions | ,|p nΨ 〉 Ψ 〉 in the field of a single regular-
ized Coulomb center whose energies may cross zero. By
making use of the wave function (130), we project the Di-
rac equation | = |h ∆Ψ〉 ε Ψ〉 on the states | pΨ 〉 and | nΨ 〉
and find the following secular equation:
det = 0,
pp pn
np nn
h h S
h S h
∆ ∆
∆ ∆
− ε − ε
− ε − ε
(132)
where = | |ijh i h j〈 〉 , | = iji j〈 〉 δ , , = ,p ni j Ψ Ψ ,
= | = |p n n pS 〈Ψ Ψ 〉 〈Ψ Ψ 〉 , | = 1〈Ψ Ψ〉 . It is easy to see
that the overlap integral
/ 2 / 2= ( ) ( ) ( ) ( )p n n p
p n
x R x RS dxdy f r g r f r g r
r r
− +
+
∫
vanishes after changing x x→ − in the first term of the bra-
ckets.
In order to calculate the coefficients ijh , it is convenient
to represent the Hamiltonian in the form
2 2 2 2
0 0= / = /p n n ph h r r h r r+ ζ + − ζ +
.
Then we obtain
0 02 2
0
2 2
0
= | | = = ,
= | | = = .
pp p p nn
n
pn p n np
p
h C h
r r
h A h
r r
ζ
ε + 〈Ψ Ψ 〉 ε + ζ −
+
ζ
−〈Ψ Ψ 〉 −ζ
+
(133)
We find that the Coulomb integral C equals
2 2 2 2
00 0
1 4= | | =
( )n
rdrC p p
r r r R r
∞
〈 〉 ×
+ + +
∫
( )2 2
2 2
0
4K ( ) ( ) ,
( )
rR f r g r
r R r
× +
+ +
(134)
where K( )k is the complete elliptic integral of the first
kind. We note that the Coulomb integral C is positive def-
inite and monotonously decreases with increasing R . Fur-
ther, the resonance integral A equals
2 2
0
1= | | =
p
A p n
r r
〈 〉
+
2 2 2 2
0 0 0
/ 2 1 1= 2 ( ) ( ) .p n
n p n
x Rdx dy f r g r
r r r r r
∞ ∞
−∞
− − + +
∫ ∫ (135)
The Coulomb integral C and resonance integral A can
be computed numerically with the functions f and g found
Fig. 10. (Color online) The energy of the electron bound state
with = 1/ 2j in the regularized Coulomb potential with Ze±
charges as a function of ζ for different values of the regulariza-
tion parameter: 0 = 0.05r R∆ (1), 0 = 0.01r R∆ (2), 0 = 0r (3).
Low Temperature Physics/Fizika Nizkikh Temperatur, 2018, v. 44, No. 5 515
E.V. Gorbar, V.P. Gusynin, and O.O. Sobol
for an isolated impurity problem. At small 0R we have
that = constC while 0A → . Asymptotically at large R they
behave as 1/C R and 2
0exp ( 1 )A R− − ε .
Finally, we obtain the energy levels
2 2 2 2 2
0= ( ) ( ) = ( ) ,pp pnh h C A∆ε ± + ± ε + ζ + ζ (136)
which are obviously symmetric with respect to the replace-
ment ∆ ∆ε → −ε in accord with the charge conjugation sym-
metry of the problem under consideration. We note that
the energy levels never cross in agreement with the avoided
crossing theorem [68]. Since the Coulomb and resonance
integrals, C and A, tend to zero as R →∞, the energy of the
system for large distances between the impurities tends to
0| |∆ε → ± ε as expected. The coefficients of the wave func-
tion of the negative energy level are given by
2 2 2 2
= ,
( ) ( ( ) ( ) )
pn
p
np pp pp pn
h
v
h h h h
−
+ + +
(137)
2 2
2 2 2 2
( ) ( )
= .
( ) ( ( ) ( ) )
pp pp pn
n
np pp pp pn
h h h
v
h h h h
+ +
+ + +
(138)
We plot the energy levels of the system for 0 = 0.05r R∆
in Fig. 11 as functions of R for = 0.65ζ , = 0.8ζ and
= 0.9ζ . In the first case, we have < = 0.7cζ ζ and the
bound state levels monotonously converge to each other as
R increases and never cross. For the couplings = 0.8ζ and
= 0.9ζ which are larger than = 0.7cζ , their behavior is no
longer monotonous. For small R , the levels converge like
in the previous case. However, after the maximal conver-
gence of the levels they go away with the subsequent in-
crease of R . This behavior is typical for the avoided cross-
ing [68].
Equation (136) and the facts that C and A monotonous-
ly depend on R and 0C ≥ imply that the level repulsion
can take place only for 0 < 0ε and the energy levels con-
verge most closely for 0=| |Cζ ε from which we can de-
termine the corresponding distance mR . Exactly at this
distance mR we have = = 0pp nnh h from Eq. (133), hence
= = 1/ 2p nv v and, consequently, the probability to find
the electron near the positively and negatively charged
impurities is the same. Furthermore, for 0| |Cζ ≈ ε , the dif-
ference of the coefficients pv and nv squared equals
2 2
2 2
2 2 2 2
2 ( ) ( )
=
( ) ( ( ) ( ) )
pp pp pp pn
n p
np pp pp pn
h h h h
v v
h h h h
+ +
−
+ + +
0sgn ( ) = sgn ( ).pph Cε + ζ (139)
The following qualitative picture appears for 0 < 0ε .
For small R , = 0pn nph h , hence | | | |n pv v and, there-
fore, the electron wave function of the negative energy
level in the LCAO method is localized mainly on the nega-
tively charged impurity. Although this result seems to be
counter-intuitive, it is quite natural. The point is that the
energy spectrum of the system for = 0R is composed of
the upper and lower continua. Since the chemical potential
in neutral graphene is zero, the electron states of the lower
continuum are occupied. For small R, the positively charged
impurity produces electron bound states which descend
from the upper continuum. Obviously, there are also charged
conjugated states localized near the negatively charged im-
purity, which rise from the lower continuum as R increases.
These states are occupied for sufficiently small R . Since
| | = | |n pv v at the point of maximal convergence = mR R ,
the probability to find the electron near the negatively and
positively charged impurities is then equal. As the distance
between impurities R increases further, the difference of
the square moduli (139) changes sign because the Coulomb
integral C decreases with increasing R . This means that
the electron wave function changes its localization to the
positively charged impurity. This change of the wave func-
Fig. 11. (Color online) The energy of the bound state levels as
functions of distance in the LCAO method for 0 = 0.05r R∆ and
= 0.65ζ (1), = 0.8ζ (2), and = 0.9ζ (3).
Fig. 12. The dependence of coefficients (137), (138) on the dis-
tance R between the impurities for the negative–energy branch of
the spectrum for = 0.85ζ : | |pv — decreasing line, | |nv — in-
creasing line. The point of the change of localization mR is
marked by dashed line.
516 Low Temperature Physics/Fizika Nizkikh Temperatur, 2018, v. 44, No. 5
Electron states in the field of charged impurities in two-dimensional Dirac systems
tion localization (the relocalization or “migration” effect) is
explicitly shown in Fig. 12 for = 0.85ζ when 0 = 0.453 < 0ε − .
The value of mR , in general, depends on 0r and on the value
of the coupling constant ζ . As the coupling increases, mR
decreases. A similar phenomenon of the change of locali-
zation of wave function takes place in the fission of quar-
konium resonances consisting of the heavy quark and anti-
quark [4].
It is clear that if 0 > 0ε , then nothing interesting happens.
Indeed, since > 0pph for 0 > 0ε , we find that | | > | |n pv v
for any distance R between the impurities. Therefore, the
wave function is always localized on the negatively charged
impurity and, consequently, the wave function of the highest
occupied state does not change its localization. Note that
the behavior of the energy levels in this case given by lines (1)
in Fig. 11 as a function of R is monotonous unlike the case
where 0ε is negative.
In Ref. 130 we considered the numerical Galerkin–Kan-
torovich variational method and obtained qualitatively sim-
ilar results.
7. Coulomb center problem in bilayer graphene
The Bernal stacked bilayer graphene forms a very inter-
esting two-dimensional physical system whose electron ex-
citations are described by the Hamiltonian with non-trivial
chiral properties [131,132]. These chiral fermions do not have
analogues in high energy physics because their energy spec-
trum 2| |E ± p
is parabolic in the two band model, like in
non-relativistic systems described by the Schrödinger equa-
tion. Still since there are positive and negative energy bands
related through the charge conjugation, this suggests that
the instability connected with diving of the lowest energy
electron bound state into the lower band may take place in
gapped bilayer graphene too. In addition, it was shown in
Ref. 133 that the chirality of quasiparticles in gapless bi-
layer graphene results in cloaked bound states that do not
hybridize with the hole continuum.
Actually, there is more to this. The parabolic spectrum
in bilayer graphene is more soft than in monolayer graphene.
Therefore, the electron-electron interactions are effectively
enhanced and may open a gap connected with the conden-
sation of electron-hole pairs in bilayer graphene. This con-
clusion agrees with the experimental data. While no gap is
observed in monolayer graphene at the neutrality point, a
small gap 2 meV is realized in bilayer graphene [134–136]
in the absence of external electromagnetic fields. A much
larger gap (~42 meV) is observed in high-mobility ABC-
stacked trilayer graphene [137] where electron excitations
have a softer dispersion 3( ) | |E ±p p . As we mentioned
above, the supercritical instability in the Coulomb center is
a precursor of the excitonic instability and gap opening in
the electron spectrum of a many-body problem. Reversing
this argument and taking into account the fact that bilayer
graphene is gapped due to the electron-electron interac-
tions, one may expect that the supercritical instability for
the electron in the field of Coulomb center should take place
in bilayer graphene.
The following qualitative consideration is important for
the analysis of the supercritical instability in bilayer gra-
phene. Since the quasiparticle kinetic energy scales like
21/ r with distance from the charged impurity in bilayer gra-
phene and is negligible compared to the Coulomb interac-
tion 1/ r , this suggests the absence of the critical charge for
the formation of a bound state in this material. On the other
hand, since the kinetic energy at small distances is larger
than the Coulomb interaction energy, this implies that the
fall-to-center should not take place in bilayer graphene.
Therefore, a priori, one may expect that the supercritical
instability in bilayer graphene should not be necessarily
related to the phenomenon of the fall-to-center. The situa-
tion is different in monolayer graphene where both the ki-
netic energy and Coulomb interaction scale equally with
distance. In fact, this equal scaling is the physical reason
why the supercritical instability in monolayer graphene is
related to the phenomenon of the fall-to-center where the
electron wave function shrinks towards the impurity as its
charge increases.
Clearly, the above heuristic reasoning is based on the
scaling in the effective low-energy theory in bilayer gra-
phene. However, the parabolic spectrum in bilayer gra-
phene [131] is valid only up to momenta 1| | / (2 )F≈ γp v ,
where 1 = 0.39eVγ is the interlayer hopping amplitude.
For larger momenta, the energy dispersion in linear as in
the monolayer graphene. Therefore, it is possible, in prin-
ciple, that as the charge of impurity increases the wave
function of the electron bound state becomes more local-
ized in the vicinity of the impurity such that the parabolic
energy spectrum does not apply. If this happens, then the
supercritical instability in bilayer graphene would proceed
like in monolayer graphene and would be connected with
the phenomenon of the fall-to-center. Actually, this scenar-
io was advocated in Ref. 138, where it was argued by using
the semiclassical approach that the supercritical instability
is an “ultra-relativistic” effect in gapless bilayer graphene
and, therefore, the value of the critical charge is the same
as in monolayer graphene.
In order to see whether the supercritical instability is in-
deed related to the atomic collapse in bilayer graphene, the
three of us studied in Ref. 139 the electron bound states in
the field of a charged impurity in bilayer graphene by us-
ing the two continuum models: the low-energy two-band
model as well as the four-band model. The main principal
advantage of the four-band model for the study of super-
critical instability in bilayer graphene is that it takes into
account the evolution of the electron energy dispersion from
the low-energy quadratic to high energy linear energy dis-
persion. Since the trigonal warping is small in bilayer gra-
phene, it could be neglected. On the other hand, the screen-
ing effects play a very essential role in view of the finite
density of states at zero energy in bilayer graphene and
Low Temperature Physics/Fizika Nizkikh Temperatur, 2018, v. 44, No. 5 517
E.V. Gorbar, V.P. Gusynin, and O.O. Sobol
should be taken into account. It will be shown that while
the lowest energy bound state in the Coulomb center in
gapped bilayer graphene indeed dives into the lower con-
tinuum, the wave function of the electron bound state does
not shrink toward the impurity as its charge increases. This
means that the supercritical instability in bilayer graphene
is not related to the phenomenon of the fall-to-center.
7.1. Two-band model
The two-band Hamiltonian, which describes low-energy
electron excitations in gapped bilayer graphene in the field
of an impurity with charge Ze, reads [131]
22
21
2
0 ( ) 1 0
= ( ),
0 1( ) 0
e( ) = ,
F p
H V r
p
ZV r
r
−
+
+ ∆ + −γ
−
κ
v
(140)
where = x yp p ip± ± and = i−p is the two-dimensional
momentum operator, 1 0.39 eVγ ≈ is the strongest inter-
layer coupling between pairs of orbitals that lie directly
below and above each other, and κ is the dielectric constant
which we choose equal to 4 in our analysis.
The energy spectrum of the free part of Hamiltonian (140)
is given by 2 2 2*( ) = ( / 2 )E m± + ∆p p , so that we have a
gap 2∆ between the lower and upper continua (the valence
and conductivity bands), where 2* 1= / 2 0.054F em mγ ≈v is
the quasiparticle mass and em is the mass of the electron.
Although the kinetic part of Hamiltonian (140) has the same
matrix structure as in monolayer graphene leading to chi-
rality and the particle-hole symmetry characteristic of rela-
tivistic systems, the dispersion relation for quasiparticles in
bilayer graphene is quadratic like in non-relativistic sys-
tems. Thus, the low-energy Hamiltonian (140) non-trivially
combines relativistic and non-relativistic properties. Clear-
ly, while a regularized potential is crucially needed for the
study of a charge instability in monolayer graphene, it suf-
fices to use the standard potential 1/ r for the two-band
model of bilayer graphene with the parabolic dispersion of
quasiparticles where the kinetic energy dominates over the
potential one at small distance.
As we discussed in the Introduction and Sec. 2, there is
a critical charge of impurity in monolayer graphene above
which the atomic collapse occurs. It is interesting what hap-
pens in bilayer graphene in view of the quadratic disper-
sion relation in this material. Although the two-band model
is applicable only up to the energies of order 1 / 4γ when
the next band becomes important, it still makes sense to study
how the supercritical Coulomb center instability is realized
in the low-energy effective two-band model. We consider
this problem in this section and discuss briefly the results
obtained in the four-band model in the next subsection.
It is convenient to define the coupling constant as
1= ( / ) /gZξ α κ γ ∆ and express distances in terms of a
characteristic length 1= /F∆λ γ ∆v and wave vectors in
its inverse, i.e., = ∆λk k (in what follows we will omit
tilde over dimensionless momenta). The eigenstates of Ha-
miltonian (140) for = ( , )TΨ χ ϕ are determined in momen-
tum space by the following system of equations:
2
2
eff2
2
2
eff2
( ) ( ) (1 ) ( ) ( ) ( ) = 0,
(2 )
( ) ( ) (1 ) ( ) ( ) ( ) = 0,
(2 )
x y
x y
d qk ik V
d qk ik V
∆
∆
− ϕ + − ε χ + χ −
π
+ χ − + ε ϕ + ϕ − π
∫
∫
k k q k q
k k q k q
(141)
where = /E∆ε ∆ and effV describes the screened potential
of impurity with charge Ze. It is defined by an analog of
the Poisson equation of the form
(2) 2
2 eff eff
4
( ) = 2 ( ) ( ) ( ),g
DV d V
πα
−∆ − πξδ − Π −
κ ∫x x y x y y
(142)
where ( )Π x is the static polarization function in bilayer
graphene. In momentum space, we easily find the follow-
ing solution to Eq. (142):
eff
2 1( ) = 2 ( ) ,
4
( )g
V q V q
q
q q
πξ
− ≡ − πξ − δ πα + Π
κ
(143)
where
4 ( )( ) =
4
( )
g
g
qV q
q q q
πα Π
δ
πακ
+ Π κ
(144)
is the correction to the Coulomb interaction due to the screen-
ing effects. The one-loop polarization function as an integral
over momentum was derived in Ref. 140 and equals
___________________________________________________
2 2 2 2
2
1
2
2( 1)
2 2 2 2
( ) = ,
( )(2 )
d k − +
+ − + −
− + − +
ε ε − − − + − + − γ Π
∆ ε ε ε + επ∫
q q q qk k k k
q
4
= 1 .
2±ε + ±
qk (145)
518 Low Temperature Physics/Fizika Nizkikh Temperatur, 2018, v. 44, No. 5
Electron states in the field of charged impurities in two-dimensional Dirac systems
For 1q , the polarization function is proportional to 2q
and tends to zero for 0q → , therefore, ( )V qδ is not singu-
lar at = 0q . For 1q , 1( ) = / ln 2 /qΠ γ ∆ π.
In polar coordinates, the angle and the absolute value of
momentum variables can be separated. In order to do this,
we note that the total z component of the pseudospin-or-
bital momentum commutes with the Hamiltonian (140)
= = ,z z z zJ L S i ∂
+ − + σ
∂θ
ˆ[ , ] = 0.zJ H (146)
Thus, we seek the spinor function in the form
( 1)
( 1)
( ) e( )
= = ,
( ) ( ) e
i j
j
i j
j
a k
b k
− θ
+ θ
χ Ψ ϕ
k
k
(147)
where j is the total angular momentum which is integer in
bilayer graphene. Substituting the spinor (147) into the
system (141) and integrating over angle, we obtain
___________________________________________________
2
1 1
0
2
1 1
0
( ) ( ) ( )[ ( , ) ( , )] = ( ),
2
( ) ( ) ( )[ ( , ) ( , )] = ( ),
2
j j j j j j
j j j j j j
k b k a k dqqa q K k q V k q a k
k a k b k dqqb q K k q V k q b k
Λ
− −
Λ
+ +
ξ + − − δ ε
π
ξ
− − − δ ε π
∫
∫
(148)
______________________________________________
where
2 2 2
| | 1/22 2
0
cos( ) 2( , ) = = ,
22 cos
j j
d j k qK k q Q
kqkqk q kq
π
−
θ θ +
+ − θ
∫
(149)
2
2 2
0
( , ) = 2 cos cos ( ).jV k q d V k q kq j
π
δ θδ + − θ θ
∫ (150)
Here ( )Q zν is the Legendre function of the second kind.
The kernels ( , )jK k q can be expressed in terms of the full
elliptic integrals of the first and second kind. Unfortunate-
ly, it is not possible to find analytically a solution to the
above system of equations. Its numerical solutions will be
given in the next subsection. However, in Ref. 139 the var-
iational method was applied in order to have an analytic
insight into the problem. It is based on the squared Schrö-
dinger equation in coordinate space. The main result is that
the critical charge of impurity crZ for the unscreened Cou-
lomb potential tends to zero as ∆ as 0∆ → , while in the
case of the screened interaction potential it tends to a finite
value for = 0∆ .
7.2. Numerical results
In order to find a numerical solution to the system of
Eqs. (148), we split the momentum interval (0, )Λ in N
equal intervals and approximate the integrals by the sums
of values of integrands at the ends of intervals multiplied
by the weight function iw of the Newton–Cotes formula of
the fifth order
0
[ ( , ) ( , )] ( ) =j j jdqq K k q V k q b q
Λ
− δ∫
=0
= [ ( , ) ( , )] ( ).
N
i i j i j i j i
i
q w K k q V k q b q− δ∑ (151)
The kernel ( , )jK k q defined in Eq. (149) has a logarith-
mic singularity at =q k . In order to deal with this singulari-
ty, we use the following regularization for the first term in
the square brackets on the left-hand side of Eq. (151):
0
( , ) ( ) =j jdqqK k q b q
Λ
∫
=0 0
= ( , ) [ ( ) ( )] ( ) ( , ).
N
i i j i j i j j j
i
q w K k q b q b k b k dqqK k q
Λ
− +∑ ∫
(152)
The last integral in Eq. (152) is not singular and could
be expressed through the elliptic integrals and generalized
hypergeometric functions by using formulas from Sec. 5.11
in Ref. 46. For = 0, 1, 2j , and 3, the corresponding expres-
sions are listed in Ref. 139.
The quadrature step h corresponds to the inverse of the
size of a graphene disc, therefore, the energy spectrum of
the system is discrete and the upper and lower continua
existing in an infinite system appear now as the sets of
closely situated discrete levels with the distance between
them proportional to h. The energy levels of the correspond-
ing upper and lower quasicontinua drift very slowly as
the impurity charge increases. On the other hand, the bound
levels inside the band gap shift towards the lower quasi-
continuum much faster. Therefore, there exists a critical
value of the impurity charge when the lowest energy bound
state approaches the highest energy state of the lower quasi-
continuum. We define this value as the critical charge. For
different values of total angular momentum j , we find
the different values of the critical charge (see Fig. 14(a)).
The minimal value of the critical charge is obtained for
= 1j . Therefore, all numerical computations in this section
are performed for this value of the total momentum.
Low Temperature Physics/Fizika Nizkikh Temperatur, 2018, v. 44, No. 5 519
E.V. Gorbar, V.P. Gusynin, and O.O. Sobol
The dependence of the critical coupling constant on gap
∆ is plotted in Fig. 13 in the cases where the screening
effects are absent and taken into account. Clearly, the criti-
cal charge decreases in both cases as gap decreases. Ac-
cording to our numerical calculations, the critical charge in
the absence of the screening effects is very well approxi-
mated by the function fit
cr 1= 1.88 /ζ ∆ γ . Such a gap de-
pendence excellently agrees with the conclusion obtained
from the variational method in Ref. 139. If the screening
effects are taken into account, then the critical coupling
constant does not tend to zero as 0∆ → . The results of the
numerical solution to system (148) are quite similar to
those obtained in Ref. 139 by using the variational method
and can be fit by the function
fit
cr 1= 0.36 4.32 / .ζ + ∆ γ (153)
In order to quantify the localization properties of the wave
function of the electron bound state in the near-critical re-
gime and check the consistency of the use of the two-band
model for the study of the supercritical instability in bilayer
graphene, we plot for crZ Z≈ in Fig. 14(b) the square of
the wave function in the two-band model in momentum
representation 2 2( ) = 2 [ ( ) ( )] /j jW k k a k b k Nπ + multiplied
by the weight factor 2 kπ , where N is the normalization
constant 2 2
0
= 2 [ ( ) ( )]j jN k a k b k dk
Λ
π +∫ . Obviously, if the
wave function is localized in momentum space in the re-
gion of momenta k Λ , then the use of the low-energy
model for the description of the supercritical phenomena
is consistent. According to Fig. 14(b), the maximum of
the wave function corresponds to 0.075k ≈ . This value is
10 times less than the cutoff of the two-band model. This
result suggests that the low-energy two-band model con-
sistently describes the supercritical behavior in bilayer
graphene.
Although our results show that the analysis of the su-
percritical instability is consistent in the two-band model, it
is still necessary to study the supercritical instability in the
four-band model of bilayer graphene. The point is that the
quadratic energy dispersion in the two-band model is re-
placed by the linear energy dispersion as in monolayer
graphene for momenta larger than 1 / (4 )Fγ v . Therefore, it
is possible, in principle, that the supercritical instability
will be affected by the electron dynamics at short distances.
The four-band model whose energy dispersion smoothly
interpolates between the low-energy quadratic and high-
energy linear in momentum energy dispersion allows us to
investigate whether the conclusions made in the present
section survive in the four-band model.
In Ref. 139 we studied the Coulomb center problem in
bilayer graphene in the four-band model. Figure 14(b) com-
pares the localization properties of the wave function of
the electron bound state in the two- and four-band models.
Clearly, the wave functions are localized in the same domain
of momentum space in these models. This proves that the
low-energy two-band model gives fairly accurate results for
the Coulomb center problem in the case of small gaps.
Fig. 13. The critical coupling constant as a function of dimen-
sionless gap 1/∆ γ with (solid line) and without (dashed line) the
screening effects taken into account.
Fig. 14. (Color online) The critical coupling constant as a func-
tion of total angular momentum j for = 3∆ meV (a). The square
of the wave function of the electron bound state for the near-cri-
tical charge = / = 0.17gZζ α κ and the gap = 3∆ meV are plotted
for the two-band and four-band models as a function of dimen-
sionless momentum, which is expressed in terms of 11 / γλ (b).
520 Low Temperature Physics/Fizika Nizkikh Temperatur, 2018, v. 44, No. 5
Electron states in the field of charged impurities in two-dimensional Dirac systems
8. Summary
In this review, we considered the electron states in the
field of charged impurities in single layer graphene as well
as bilayer graphene paying the special attention to the phe-
nomenon of supercritical instability and its realization in
the presence of several charged impurities and the external
magnetic field. Such states are physically relevant for the
actively developing area of graphene quantum dots whose
shape and boundary conditions at edges strongly affect the
confined states [141]. In addition, the electron states in the
field of charged impurities are crucially important for the
transport properties of graphene [142].
In the supercritical Coulomb center problem, we showed
that the “fall-to-center” phenomenon arises if Zα exceeds
the critical value 1/2 leading to the appearance of quasista-
tionary levels with complex energies. The energy of quasi-
stationary states in the case of gapless quasiparticles has
a characteristic essential-singularity type dependence on
the coupling constant reflecting the scale invariance of the
Coulomb potential. It is found that a quasiparticle gap sta-
bilizes the system decreasing the imaginary part | Im |E of
quasistationary states.
We reviewed also the experiments in which the super-
critical behavior was observed. The STM measurements
revealed the formation of resonances around artificial nuclei
(clusters of charged calcium dimers) fabricated on gated gra-
phene devices via atomic manipulation techniques. In a mag-
netic field, the STM and Landau level spectroscopy measu-
rements demonstrate that the strength of the impurity can
be tuned by controlling the occupation of Landau-level
states with a gate-voltage. At low occupation the impurity
is screened becoming essentially invisible. The screening
of the potential of the impurity diminishes as states are filled
until, for fully occupied Landau levels, the unscreened im-
purity significantly perturbs the spectrum in its vicinity. In
this regime it is possible to observe the Landau-level split-
ting into discrete states due to lifting the orbital degeneracy.
Further, we studied theoretically the electron states in
the field of a charged impurity in graphene in a magnetic
field. The charged impurity removes the degeneracy of Lan-
dau levels converting them into band like structures. As the
charge of impurity grows, the repulsion of sublevels of dif-
ferent Landau levels with the same value of orbital momen-
tum takes place leading to the redistribution of the wave
function profiles of these sublevels near the impurity [62,65].
This qualitatively corresponds to the formation of reso-
nance states in the traditional version of supercritical insta-
bility. Taking into account the polarization effects in a mag-
netic field allowed us to explain the tuning of the effective
charge of the impurity by the gate voltage in agreement
with the recent experiments.
In addition, the two-particle bound state problem for
gapped graphene in the presence of a Coulomb impurity
was studied. A variational approach, using the projected Ha-
miltonian and Chandrasekhar–Dirac spinors as trial wave
functions, predicts the existence of at least one bound state.
In contrast to the Schrödinger case, the variational energy
functional is not a homogeneous function of the coupling
constant α. As a consequence, the optimal values of the va-
riational parameters depend on α and the optimal binding
energy has a more complicated functional dependence on α.
In particular, the binding energy increases with respect to
the nonrelativistic case.
We studied also the supercritical instability in gapped gra-
phene with two charged impurities separated by distance R
in the case where the charges of impurities are subcritical,
whereas their total charge exceeds a critical one. The criti-
cal distance crR in the system of two charged centers is
defined as that at which the electron bound state with the
lowest energy reaches the boundary of the lower continu-
um. Since the variables in the Dirac problem with two
Coulomb centers are not separable in any known orthogo-
nal coordinate system, this problem does not admit an ana-
lytic solution. Therefore, the critical distance crR could be
calculated only with the help of approximate methods. The
LCAO technique and variational Kantorovich method give
quite similar results. They show that the critical distance
crR increases as the quasiparticle gap decreases. The transi-
tion to the supercritical regime is signaled by the appear-
ance of quasistationary states in the lower continuum.
A new type of supercritical behavior in gapped gra-
phene with two oppositely charged impurities was revealed
by studying the two-dimensional Dirac equation for quasi-
particles with the Coulomb potential regularized at small
distances. By utilizing the technique of linear combination
of atomic orbitals and the variational Galerkin–Kantoro-
vich method, it was shown that for supercritical electric
dipole the wave function of the electron bound state changes
its localization from the negatively charged impurity to the
positively charged one as the distance between the impuri-
ties changes. Such a migration of the wave function corre-
sponds to the electron and hole spontaneously created from
the vacuum in bound states screening the positively and
negatively charged impurities of the supercritical electric
dipole, respectively. The obtained results were generalized
in Ref. 130 to a particle-hole asymmetric case, where the
charges of impurities differ in signs and absolute values,
and it was demonstrated that the necessary energetic condi-
tion for the supercriticality of novel type to occur is that
the energy levels of single positively and negatively charg-
ed impurities traverse together the energy distance separat-
ing the upper and lower continua.
Finally, the supercritical instability in gapped bilayer
graphene was studied in the low energy two-band as well
as four-band continuum models. The different scalings of
the kinetic energy of quasiparticles and the Coulomb inter-
action with respect to the distance to the charged impurity
ensure that the wave function of the electron bound state
does not shrink toward the impurity as its charge increases.
Low Temperature Physics/Fizika Nizkikh Temperatur, 2018, v. 44, No. 5 521
E.V. Gorbar, V.P. Gusynin, and O.O. Sobol
This results in the absence of the fall-to-center phenome-
non in bilayer graphene although the supercritical instabil-
ity is realized. It was found that the screening effects are
crucially important in bilayer graphene. If they are neglect-
ed, then the critical value for the impurity charge as the
lowest energy bound state dives into the lower continuum
tends to zero as the gap ∆ vanishes. If the screened Coulomb
interaction is considered, then the critical charge tends to
a finite value when gap goes to zero.
The work of E.V.G. and V.P.G. is partially supported
by the National Academy of Sciences of Ukraine (project
0116U003191) and by its Program of Fundamental Re-
search of the Department of Physics and Astronomy (pro-
ject No. 0117U000240). V.P.G. acknowledges the support
of the RISE Project CoExAN GA644076. References
_______
1. L. Landau and E. Lifshitz, Quantum Mechanics. Non-relati-
vistic Theory, Butterworth-Heinemann, Oxford (2004).
2. I.Y. Pomeranchuk and Y.A. Smorodinsky, J. Phys. USSR 9,
97 (1945).
3. Y.B. Zeldovich and V.S. Popov, Sov. Phys. Usp. 14, 673 (1972).
4. W. Greiner, B. Müller, and J. Rafelski, Quantum Electro-
dynamics of Strong Fields. With an Introduction into Mo-
dern Relativistic Quantum Mechanics, Texts and Mono-
graphs in Physics, Springer-Verlag, Berlin-Heidelberg (1985).
5. S.S. Gershtein and Y.B. Zeldovich, Sov. Phys. JETP 30, 358
(1970).
6. J. Rafelski, L.P. Fulcher, and W. Greiner, Phys. Rev. Lett.
27, 958 (1971).
7. B. Müller, H. Peitz, J. Rafelski, and W. Greiner, Phys. Rev.
Lett. 28, 1235 (1972).
8. L.A. Fal’kovsky, Sov. Phys. Usp. 11, 1 (1968).
9. V.S. Edel’man, Adv. Phys. 25, 555 (1976).
10. C. Herring, Phys. Rev. 52, 365 (1937).
11. N.P. Armitage, E.J. Mele, and A. Vishwanath, Rev. Mod.
Phys. 90, 015001 (2018).
12. A.V. Shytov, M.I. Katsnelson, and L.S. Levitov, Phys. Rev.
Lett. 99, 246802 (2007).
13. A.V. Shytov, M.I. Katsnelson, and L.S. Levitov, Phys. Rev.
Lett. 99, 236801 (2007).
14. V.M. Pereira, J. Nilsson, and A.H. Castro Neto, Phys. Rev.
Lett. 99, 166802 (2007).
15. M.M. Fogler, D.S. Novikov, and B.I. Shklovskii, Phys. Rev.
B 76, 233402 (2007).
16. Y. Wang, D. Wong, A.V. Shytov, V.W. Brar, S. Choi, Q. Wu,
H.-Z. Tsai, W. Regan, A. Zettl, R.K. Kawakami, S.G. Louie,
L.S. Levitov, and M.F. Crommie, Science 340, 734 (2013).
17. A. Feher, I.A. Gospodarev, V.I. Grishaev, K.V. Kravchenko,
E.V. Manzhelii, E.S. Syrkin, and S.B. Feodos’ev, Fiz. Nizk.
Temp. 35, 862 (2009) [Low Temp. Phys. 35, 679 (2009)].
18. V.M. Loktev and Y.G. Pogorelov, Fiz. Nizk. Temp. 38, 993
(2012) [Low Temp. Phys. 38, 792 (2012)].
19. Y.V. Skrypnyk and V.M. Loktev, Fiz. Nizk. Temp. 33, 1002
(2007) [Low Temp. Phys. 33, 762 (2007)].
20. B. Dóra, Fiz. Nizk. Temp. 34, 1020 (2008) [Low Temp. Phys.
34, 801 (2008)].
21. O.V. Gamayun, E.V. Gorbar, and V.P. Gusynin, Phys. Rev.
B 80, 165429 (2009).
22. J. Wang, H.A. Fertig, and G. Murthy, Phys. Rev. Lett. 104,
186401 (2010).
23. J. Sabio, F. Sols, and F. Guinea, Phys. Rev. B 81, 45428
(2010).
24. D.V. Khveshchenko, Phys. Rev. Lett. 87, 246802 (2001).
25. E.V. Gorbar, V.P. Gusynin, V.A. Miransky, and I.A. Shovkovy,
Phys. Rev. B 66, 045108 (2002).
26. E.V. Gorbar, V.P. Gusynin, V.A. Miransky, and I.A. Shovkovy,
Phys. Lett. A 313, 472 (2003).
27. D.V. Khveshchenko and H. Leal, Nucl. Phys. B 687, 323
(2004).
28. O.V. Gamayun, E.V. Gorbar, and V.P. Gusynin, Phys. Rev. B
81, 75429 (2010).
29. J. González, Phys. Rev. B 85, 85420 (2012).
30. J.E. Drut and T.A. Lähde, Phys. Rev. Lett. 102, 026802 (2009).
31. J.E. Drut and T.A. Lähde, Phys. Rev. B 79, 241405 (2009).
32. W. Armour, S. Hands, and C. Strouthos, Phys. Rev. B 81,
125105 (2010).
33. P.V. Buividovich and M.I. Polikarpov, Phys. Rev. B 86,
245117 (2012).
34. P.I. Fomin, V.P. Gusynin, V.A. Miransky, and Y.A. Sitenko,
Riv. Nuovo Cimento 6, 1 (1983).
35. R.D. Peccei, Nature 332, 492 (1988).
36. A. Luican-Mayer, M. Kharitonov, G. Li, C.-P. Lu, I. Skachko,
A.-M.B. Gonçalves, K. Watanabe, T. Taniguchi, and E.Y.
Andrei, Phys. Rev. Lett. 112, 036804 (2014).
37. J. Mao, Y. Jiang, D. Moldovan, G. Li, K. Watanabe, T.
Taniguchi, M.R. Masir, F.M. Peeters, and E.Y. Andrei,
Nature Phys. 12, 545 (2016).
38. G. Giovannetti, P.A. Khomyakov, G. Brocks, P.J. Kelly, and
J. van den Brink, Phys. Rev. B 76, 073103 (2007).
39. Y.-W. Son, M.L. Cohen, and S.G. Louie, Phys. Rev. Lett. 97,
216803 (2006).
40. D.S. Novikov, Phys. Rev. B 76, 245435 (2007).
41. V.M. Pereira, V.N. Kotov, and A.H. Castro Neto, Phys. Rev.
B 78, 085101 (2008).
42. I.S. Terekhov, A.I. Milstein, V.N. Kotov, and O.P. Sushkov,
Phys. Rev. Lett. 100, 076803 (2008).
43. A. Shytov, M. Rudner, N. Gu, M. Katsnelson, and L. Levitov,
Solid State Commun. 149, 1087 (2009).
44. A.H. Castro Neto, V.N. Kotov, J. Nilsson, V.M. Pereira,
N.M.R. Peres, and B. Uchoa, Solid State Commun. 149,
1094 (2009).
45. C.G. Darwin, Philos. Mag. 25, 201 (1913).
46. I.S. Gradshtein and I.M. Ryzhik, Table of Integrals, Series,
and Products, Academic, Orlando, FL (1980).
47. V.R. Khalilov and C.-L. Ho, Mod. Phys. Lett. A 13, 615 (1998).
48. P. Fomin, V. Gusynin, and V. Miransky, Phys. Lett. B 78,
136 (1978).
49. V.S. Popov, Sov. Phys. JETP 32, 526 (1971).
50. P. Fomin and V. Miransky, Phys. Lett. B 64, 166 (1976).
522 Low Temperature Physics/Fizika Nizkikh Temperatur, 2018, v. 44, No. 5
https://doi.org/10.1070/PU1972v014n06ABEH004735
https://doi.org/10.1103/PhysRevLett.27.958
https://doi.org/10.1103/PhysRevLett.28.1235
https://doi.org/10.1103/PhysRevLett.28.1235
https://doi.org/10.1070/PU1968v011n01ABEH003721
https://doi.org/10.1080/00018737600101452
https://doi.org/10.1103/PhysRev.52.365
https://doi.org/10.1103/PhysRevLett.99.246802
https://doi.org/10.1103/PhysRevLett.99.246802
https://doi.org/10.1103/PhysRevLett.99.236801
https://doi.org/10.1103/PhysRevLett.99.236801
https://doi.org/10.1103/PhysRevLett.99.166802
https://doi.org/10.1103/PhysRevLett.99.166802
https://doi.org/10.1103/PhysRevB.76.233402
https://doi.org/10.1103/PhysRevB.76.233402
https://doi.org/10.1126/science.1234320
https://doi.org/10.1063/1.3224726
https://doi.org/10.1063/1.4743562
https://doi.org/10.1063/1.2780170
https://doi.org/10.1063/1.2981390
https://doi.org/10.1103/PhysRevB.80.165429
https://doi.org/10.1103/PhysRevB.80.165429
https://doi.org/10.1103/PhysRevLett.104.186401
https://doi.org/10.1103/PhysRevB.81.045428
https://doi.org/10.1103/PhysRevLett.87.246802
https://doi.org/10.1103/PhysRevB.66.045108
https://doi.org/10.1016/S0375-9601(03)00846-6
https://doi.org/10.1016/j.nuclphysb.2004.03.020
https://doi.org/10.1103/PhysRevB.81.075429
https://doi.org/10.1103/PhysRevB.85.085420
https://doi.org/10.1103/PhysRevLett.102.026802
https://doi.org/10.1103/PhysRevB.79.241405
https://doi.org/10.1103/PhysRevB.81.125105
https://doi.org/10.1103/PhysRevB.86.245117
https://doi.org/10.1007/BF02740014
https://doi.org/10.1038/332492b0
https://doi.org/10.1103/PhysRevLett.112.036804
https://doi.org/10.1103/PhysRevB.76.073103
https://doi.org/10.1103/PhysRevLett.97.216803
https://doi.org/10.1103/PhysRevB.76.245435
https://doi.org/10.1103/PhysRevB.78.085101
https://doi.org/10.1103/PhysRevB.78.085101
https://doi.org/10.1103/PhysRevLett.100.076803
https://doi.org/10.1016/j.ssc.2009.02.043
https://doi.org/10.1016/j.ssc.2009.02.040
https://doi.org/10.1142/S0217732398000668
https://doi.org/10.1016/0370-2693(78)90366-0
https://doi.org/10.1016/0370-2693(76)90321-X
Electron states in the field of charged impurities in two-dimensional Dirac systems
51. Y. Wang, V.W. Brar, A.V. Shytov, Q. Wu, W. Regan, H.-Z.
Tsai, A. Zettl, L.S. Levitov, and M.F. Crommie, Nature
Phys. 8, 653 (2012).
52. J. Li, W.-D. Schneider, and R. Berndt, Phys. Rev. B 56, 7656
(1997).
53. C. Wittneven, R. Dombrowski, M. Morgenstern, and R.
Wiesendanger, Phys. Rev. Lett. 81, 5616 (1998).
54. I.V. Krive and S.A. Naftulin, Phys. Rev. D 46, 2737 (1992).
55. K.G. Klimenko, Theor. Math. Phys. 89, 1161 (1991).
56. V.P. Gusynin, V.A. Miransky, and I.A. Shovkovy, Phys.
Rev. Lett. 73, 3499 (1994).
57. V. Gusynin, V. Miransky, and I. Shovkovy, Nucl. Phys. B
462, 249 (1996).
58. E.V. Gorbar, V.P. Gusynin, V.A. Miransky, and I.A. Shovkovy,
Phys. Scripta T146, 14018 (2012).
59. V.N. Oraevskii, A.I. Rex, and V.B. Semikoz, Sov. Phys.
JETP 45, 428 (1977).
60. P. Schluter, G. Soff, K.H. Wietschorke, and W. Greiner, J.
Phys. B: At. Molec. Phys. 18, 1685 (1985).
61. C.-L. Ho and V.R. Khalilov, Phys. Rev. A 61, 32104 (2000).
62. O.V. Gamayun, E.V. Gorbar, and V.P. Gusynin, Phys. Rev. B
83, 235104 (2011).
63. T. Maier and H. Siedentop, J. Math. Phys. 53, 095207 (2012).
64. S.C. Kim and S.-R. Eric Yang, Ann. Phys. 347, 21 (2014).
65. D. Moldovan, M.R. Masir, and F.M. Peeters, 2D Materials
5, 015017 (2018).
66. G.W. Semenoff, I.A. Shovkovy, and L.C.R. Wijewardhana,
Phys. Rev. D 60, 105024 (1999).
67. O.O. Sobol, P.K. Pyatkovskiy, E.V. Gorbar, and V.P. Gusynin,
Phys. Rev. B 94, 115409 (2016).
68. J. von Neumann and E.P. Wigner, Phys. Z. 30, 467 (1929).
69. B. Müller, J. Rafelski, and W. Greiner, Z. Phys. 257, 62 (1972).
70. A. Kovner and B. Rosenstein, Phys. Rev. B 42, 4748 (1990).
71. N. Dorey and N. Mavromatos, Nucl. Phys. B 386, 614 (1992).
72. E.C. Marino, Nucl. Phys. B 408, 551 (1993).
73. E.V. Gorbar, V.P. Gusynin, and V.A. Miransky, Phys. Rev. D
64, 105028 (2001).
74. P.K. Pyatkovskiy and V.P. Gusynin, Phys. Rev. B 83,
075422 (2011).
75. R. Roldán, J.-N. Fuchs, and M.O. Goerbig, Phys. Rev. B 80,
085408 (2009).
76. R. Roldán, M.O. Goerbig, and J.-N. Fuchs, Semicond. Sci.
Technol. 25, 034005 (2010).
77. H.A. Mizes and J.S. Foster, Science 244, 559 (1989).
78. Y. Zhang, Y. Barlas, and K. Yang, Phys. Rev. B 85, 165423
(2012).
79. A. De Martino and R. Egger, Phys. Rev. B 95, 085418 (2017).
80. J. Lee, D. Wong, J. Velasco, Jr., J.F. Rodriguez-Nieva, S. Kahn,
H.-Z. Tsai, T. Taniguchi, K. Watanabe, A. Zettl, F. Wang, L.S.
Levitov, and M.F. Crommie, Nature Phys. 12, 1032 (2016).
81. C. Gutiérrez, L. Brown, C.-J. Kim, J. Park, and A.N. Pasupathy,
Nature Phys. 12, 1069 (2016).
82. S. Chandrasekhar, Astrophys. J. 100, 176 (1944).
83. H.A. Bethe and E.E. Salpeter, Quantum Mechanics of One-
and Two-Electron Atoms, Academic Press, New York (1957).
84. R.N. Hill, J. Math. Phys. 18, 2316 (1977).
85. B.H. Bransden and C.J. Joachain, Physics of Atoms and
Molecules, Longmont Group Limited, Hong Kong (1983).
86. T. Andersen, Phys. Rep. 394, 157 (2004).
87. H. Høgaasen, J.-M. Richard, and P. Sorba, Am. J. Phys. 78,
86 (2010).
88. D.E. Phelps and K.K. Bajaj, Phys. Rev. B 27, 4883 (1983).
89. T. Pang and S.G. Louie, Phys. Rev. Lett. 65, 1635 (1990).
90. D.M. Larsen and S.Y. McCann, Phys. Rev. B 45, 3485 (1992).
91. D.M. Larsen and S.Y. McCann, Phys. Rev. B 46, 3966 (1992).
92. N.P. Sandler and C.R. Proetto, Phys. Rev. B 46, 7707 (1992).
93. M.V. Ivanov and P. Schmelcher, Phys. Rev. B 65, 205313
(2002).
94. S. Huant, S.P. Najda, and B. Etienne, Phys. Rev. Lett. 65,
1486 (1990).
95. A.J. Shields, M. Pepper, M.Y. Simmons, and D.A. Ritchie,
Phys. Rev. B 52, 7841 (1995).
96. G.E. Brown and D.G. Ravenhall, Proc. Royal Soc. A 208,
552 (1951).
97. A. Kolakowska, J.D. Talman, and K. Aashamar, Phys. Rev. A
53, 168 (1996).
98. H. Nakatsuji and H. Nakashima, Phys. Rev. Lett. 95,
050407 (2005).
99. J. Sucher, Phys. Rev. A 22, 348 (1980).
100. J. Sucher, Int. J. Quantum Chem. 25, 3 (1984).
101. W. Häusler and R. Egger, Phys. Rev. B 80, 161402 (2009).
102. R. Egger, A. De Martino, H. Siedentop, and E. Stockmeyer,
J. Phys. A: Math. Theor. 43, 215202 (2010).
103. M.S. Marinov and V.S. Popov, Sov. Phys. JETP 41, 205
(1975).
104. C. Cohen-Tannoudji, B. Diu, and F. Laloe, Quantum Me-
chanics, Vol. 2, Hermann, Paris (1977).
105. D. Klöpfer, A. De Martino, D. Matrasulov, and R. Egger,
Eur. Phys. J. B 87, 187 (2014).
106. D. Klöpfer, Critical Effects and Universality in Graphene
Quantum Dots, PhD thesis, Düsseldorf (2014).
107. V.I. Matveev, D.U. Matrasulov, and H.Y. Rakhimov, Phys.
At. Nucl. 63, 318 (2000).
108. V.V. Bondarchuk, I.M. Shvab, D.I. Bondar, and A.V.
Katernoga, Phys. Rev. A 76, 062507 (2007).
109. V.S. Popov, Sov. J. Nucl. Phys. 14, 257 (1972).
110. V.R. Khalilov, Theor. Math. Phys. 175, 637 (2013).
111. B. Chakraborty, K.S. Gupta, and S. Sen, J. Phys. A: Math.
Theor. 46, 055303 (2013).
112. B. Chakraborty, K.S. Gupta, and S. Sen, J. Phys. Conf.
Series 442, 012017 (2013).
113. V.S. Popov, JETP Lett. 16, 251 (1972).
114. V.S. Popov, Phys. At. Nucl. 64, 367 (2001).
115. M.S. Marinov, V.S. Popov, and V.S. Stolin, J. Comput.
Phys. 19, 241 (1975).
116. O.O. Sobol, Ukr. J. Phys. 59, 531 (2014).
117. O.O. Sobol, E.V. Gorbar, and V.P. Gusynin, Phys. Rev. B
88, 205116 (2013).
118. A. De Martino, D. Klöpfer, D. Matrasulov, and R. Egger,
Phys. Rev. Lett. 112, 186603 (2014).
119. D.U. Matrasulov, V.I. Matveev, and M.M. Musakhanov,
Phys. Rev. A 60, 4140 (1999).
Low Temperature Physics/Fizika Nizkikh Temperatur, 2018, v. 44, No. 5 523
https://doi.org/10.1103/PhysRevB.56.7656
https://doi.org/10.1103/PhysRevLett.81.5616
https://doi.org/10.1103/PhysRevD.46.2737
https://doi.org/10.1007/BF01015908
https://doi.org/10.1103/PhysRevLett.73.3499
https://doi.org/10.1103/PhysRevLett.73.3499
https://doi.org/10.1016/0550-3213(96)00021-1
https://doi.org/10.1088/0031-8949/2012/T146/014018
https://doi.org/10.1088/0022-3700/18/9/005
https://doi.org/10.1088/0022-3700/18/9/005
https://doi.org/10.1103/PhysRevA.61.032104
https://doi.org/10.1103/PhysRevB.83.235104
https://doi.org/10.1063/1.4728982
https://doi.org/10.1016/j.aop.2014.04.022
https://doi.org/10.1103/PhysRevD.60.105024
https://doi.org/10.1103/PhysRevB.94.115409
https://doi.org/10.1007/BF01398198
https://doi.org/10.1103/PhysRevB.42.4748
https://doi.org/10.1016/0550-3213(92)90632-L
https://doi.org/10.1016/0550-3213(93)90379-4
https://doi.org/10.1103/PhysRevD.64.105028
https://doi.org/10.1103/PhysRevB.83.075422
https://doi.org/10.1103/PhysRevB.80.085408
https://doi.org/10.1088/0268-1242/25/3/034005
https://doi.org/10.1088/0268-1242/25/3/034005
https://doi.org/10.1126/science.244.4904.559
https://doi.org/10.1103/PhysRevB.85.165423
https://doi.org/10.1103/PhysRevB.95.085418
https://doi.org/10.1086/144654
https://doi.org/10.1007/978-3-662-12869-5
https://doi.org/10.1007/978-3-662-12869-5
https://doi.org/10.1063/1.523241
https://doi.org/10.1016/j.physrep.2004.01.001
https://doi.org/10.1119/1.3236392
https://doi.org/10.1103/PhysRevB.27.4883
https://doi.org/10.1103/PhysRevLett.65.1635
https://doi.org/10.1103/PhysRevB.45.3485
https://doi.org/10.1103/PhysRevB.46.3966
https://doi.org/10.1103/PhysRevB.46.7707
https://doi.org/10.1103/PhysRevB.65.205313
https://doi.org/10.1103/PhysRevLett.65.1486
https://doi.org/10.1103/PhysRevB.52.7841
https://doi.org/10.1098/rspa.1951.0181
https://doi.org/10.1103/PhysRevA.53.168
https://doi.org/10.1103/PhysRevLett.95.050407
https://doi.org/10.1103/PhysRevA.22.348
https://doi.org/10.1002/qua.560250103
https://doi.org/10.1103/PhysRevB.80.161402
https://doi.org/10.1088/1751-8113/43/21/215202
https://doi.org/10.1140/epjb/e2014-50414-8
https://doi.org/10.1134/1.855637
https://doi.org/10.1134/1.855637
https://doi.org/10.1103/PhysRevA.76.062507
https://doi.org/10.1007/s11232-013-0052-y
https://doi.org/10.1088/1751-8113/46/5/055303
https://doi.org/10.1088/1751-8113/46/5/055303
https://doi.org/10.1088/1742-6596/442/1/012017
https://doi.org/10.1088/1742-6596/442/1/012017
https://doi.org/10.1134/1.1358463
https://doi.org/10.1016/0021-9991(75)90075-3
https://doi.org/10.1016/0021-9991(75)90075-3
https://doi.org/10.15407/ujpe59.05.0531
https://doi.org/10.1103/PhysRevB.88.205116
https://doi.org/10.1103/PhysRevLett.112.186603
https://doi.org/10.1103/PhysRevA.60.4140
E.V. Gorbar, V.P. Gusynin, and O.O. Sobol
120. K. Connolly and D.J. Griffiths, Am. J. Phys. 75, 524 (2007).
121. J.E. Turner, Am. J. Phys. 45, 758 (1977).
122. E.V. Gorbar, V.P. Gusynin, and O.O. Sobol, Europhys.
Lett. 111, 37003 (2015).
123. D.I. Abramov and I.V. Komarov, Theor. Math. Phys. 13,
1090 (1972).
124. H.E. Camblong, L.N. Epele, H. Fanchiotti, and C.A. Garca,
Canal Phys. Rev. Lett. 87, 220402 (2001).
125. D. Schumayer, B.P. van Zyl, R.K. Bhaduri, and D.A.W.
Hutchinson, Europhys. Lett. 89, 13001 (2010).
126. V. Efimov, Phys. Lett. B 33, 563 (1970).
127. E. Braaten and H.-W. Hammer, Phys. Rep. 428, 259 (2006).
128. A.O. Gogolin, C. Mora, and R. Egger, Phys. Rev. Lett. 100,
140404 (2008).
129. M. Abramowitz and I.A. Stegun, Handbook of Mathema-
tical Functions, Dover, New York (1964).
130. E.V. Gorbar, V.P. Gusynin, and O.O. Sobol, Phys. Rev. B
92, 235417 (2015).
131. E. McCann and V.I. Fal’ko, Phys. Rev. Lett. 96, 086805
(2006).
132. K.S. Novoselov, E. McCann, S.V. Morozov, V.I. Fal’ko,
M.I. Katsnelson, U. Zeitler, D. Jiang, F. Schedin, and A.K.
Geim, Nature Phys. 2, 177 (2006).
133. A.V. Shytov, Cloaked Resonant States in Bilayer
Graphene, ArXiv [cond-mat.mes-hall], 1506.02839 (2015).
134. J. Martin, B.E. Feldman, R.T. Weitz, M.T. Allen, and
A. Yacoby, Phys. Rev. Lett. 105, 256806 (2010).
135. R.T. Weitz, M.T. Allen, B.E. Feldman, J. Martin, and
A. Yacoby, Science 330, 812 (2010).
136. F. Freitag, J. Trbovic, M. Weiss, and C. Schönenberger,
Phys. Rev. Lett. 108, 076602 (2012).
137. Y. Lee, D. Tran, K. Myhro, J. Velasco, N. Gillgren, C.N.
Lau, Y. Barlas, J.M. Poumirol, D. Smirnov, and F. Guinea,
Nature Commun. 5, 5656 (2014).
138. E.B. Kolomeisky, J.P. Straley, and D.L. Abrams, J. Phys.:
Cond. Mat. 28, 47LT 01 (2016).
139. D.O. Oriekhov, O.O. Sobol, E.V. Gorbar, and V.P. Gusynin,
Phys. Rev. B 96, 165403 (2017).
140. R. Nandkishore and L. Levitov, Phys. Rev. Lett. 104,
156803 (2010).
141. R. Van Pottelberge, M. Zarenia, P. Vasilopoulos, and F.M.
Peeters, Phys. Rev. B 95, 245410 (2017).
142. S. Das Sarma, S. Adam, E.H. Hwang, and E. Rossi, Rev.
Mod. Phys. 83, 407 (2011).
___________________________
524 Low Temperature Physics/Fizika Nizkikh Temperatur, 2018, v. 44, No. 5
https://doi.org/10.1119/1.2710485
https://doi.org/10.1119/1.10767
https://doi.org/10.1209/0295-5075/111/37003
https://doi.org/10.1209/0295-5075/111/37003
https://doi.org/10.1007/BF01035530
https://doi.org/10.1103/PhysRevLett.87.220402
https://doi.org/10.1209/0295-5075/89/13001
https://doi.org/10.1016/0370-2693(70)90349-7
https://doi.org/10.1016/j.physrep.2006.03.001
https://doi.org/10.1103/PhysRevLett.100.140404
https://doi.org/10.1103/PhysRevB.92.235417
https://doi.org/10.1103/PhysRevLett.96.086805
https://doi.org/10.1103/PhysRevLett.105.256806
https://doi.org/10.1126/science.1194988
https://doi.org/10.1103/PhysRevLett.108.076602
https://doi.org/10.1038/ncomms6656
https://doi.org/10.1103/PhysRevB.96.165403
https://doi.org/10.1103/PhysRevLett.104.156803
https://doi.org/10.1103/PhysRevB.95.245410
https://doi.org/10.1103/RevModPhys.83.407
https://doi.org/10.1103/RevModPhys.83.407
1. Introduction
2. Atomic collapse in monolayer graphene
2.1. Resonance states in gapless graphene in quasiclassical approach
2.2. Supercritical instability in graphene
2.2.1. Gapped graphene, subcritical regime
2.2.2. Gapped graphene, supercritical regime
2.2.3. Gapless graphene
2.3. Experimental observation of the atomic collapse in graphene
3. Supercritical instability in a magnetic field
3.1. The Coulomb center in a magnetic field
3.2. Tuning the screening of charged impurity with chemical potential
3.3. Screening charged impurities and lifting the orbital degeneracy in graphene by populating Landau levels
4. Two-electron bound states near a Coulomb impurity
5. Two-center problem
5.1. LCAO approach for symmetric two-center problem
5.2. Variational method
5.3. Quasistationary states
6. The dipole problem and migration of the wave function
6.1. Point electric dipole in graphene
6.2. Migration of the wave function in the finite dipole potential
7. Coulomb center problem in bilayer graphene
7.1. Two-band model
7.2. Numerical results
8. Summary
References
|