### Electron states in the field of charged impurities in two-dimensional Dirac systems

### (Review Article)

### E.V. Gorbar

^{1,2}

### , V.P. Gusynin

^{2}

### , and O.O. Sobol

^{1}

*1**Department of Physics, Taras Shevchenko National University of Kiev, Kiev 03680, Ukraine *
E-mail: [email protected]

*2**Bogolyubov 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

**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 **p**^{2}/ (2 )*m** _{e}*

^{2}/ (2

*m r*

_{e}^{2}) as

*r*→0. 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/*V r* *r** ^{n}* with

*n*≥2. In fact, the Schrö- dinger equation with the potential

*V r*( ) =−β/

*r*

^{2}provides the canonical textbook example of the fall-to-center in quan- tum mechanics [1], which takes place for β>

^{2}/ (8 )

*m*

*when the energy spectrum is not bounded from below. If the interaction potential is regularized at some distance*

_{e}*r*

_{0}, then the electron wave function of the ground state is local- ized in the region of the radius

*r*

_{0}which shrinks to the origin as

*r*

_{0}→0

_{. }

Still physically as | |**p** attains the value of order *m c** _{e}* ,
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

*r*→0. 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

*Z*170 [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 α* _{g}* = e /(

^{2}

*v*

*) 2.2≈ , where*

_{F}*v*

*≈*

_{F}*c*/300 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 *V r** _{C}*( ) = e /

^{2}

*r. The crucial advantage of graphene com-*pared to QED is its experimental accessibility where atom-

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 α α> * _{c}*1 (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

*Z*

*α in the prob- lem of the Coulomb center. In the strong coupling regime*

_{c}> _{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}* _{z}* ( ),

*H* **p** ξ *v* σ**p**+ ξ∆σ +*V* **r** _{(1) }
where *v** _{F}* is the Fermi velocity of graphene,

**p**=−

*i*∇ is the canonical momentum, σ

*are the Pauli matrices, ∆ is a quasiparticle gap, and ξ is an index, which corresponds to the valley*

_{i}*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 Θ=

*is*

_{2 1}σ

*K*

_{: }

( , = 1) 1= ( , = 1),

*H* ^{−} *H*

Θ **p** ξ + Θ − ξ −**p** (2)

where *s*_{2} 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 α* _{c}* = 1/ 2. 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 *M M*< * _{c}* = e /

*Z*

^{2}

*c.*

These states can be constructed quasiclassically from relativistic dynamics described by the Hamiltonian

= * _{F}* | | ( )

*H* *v* **p** +*V r* , where *V r*( ) =−*Z*e / ( )^{2} κ*r* and κ is a di-
electric constant. The collapsing trajectories with angular
momenta *M M*< * _{c}* = e / (

*Z*

^{2}κ

*v*

*) are separated from non- falling trajectories by a centrifugal barrier. This is mani- fested in the expression for the radial momentum square*

_{F}2 2 2

2 2

2

= e .

*r* *F* *Z* *M*

*p* *E*

*r* *r*

−

+ −

κ

*v* (3)

Clearly, there is a classically forbidden region, the an-
nulus *r*_{1}< <*r r*_{2}, *r*_{1,2}= ( e /*Z* ^{2} κ*Mv** _{F}*)/ | |

*E*, 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 *r*_{0} is a regularization parameter, which is of order of
lattice spacing. Evaluating the integral with logarithmic
accuracy, we obtain γln ( e / (*Z* ^{2} κ*r E*_{0}| |)) =π*n*, where

2 2 1/2

(*M*_{c}*M* )

γ ≡ − , which gives the quasi-Rydberg states

2 /

e e0 * ^{n}* , > 0.

*n* *Z*

*E* *n*

*r*

−π γ

≈ −κ ^{} (4)

The energies of these states converge to zero, *E** _{n}* →0

_{, }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

*p*

*and obtain the tunnel- ing action*

_{r}### ( )

2 2 2

2 1

= = .

*r* *c*

*F* *c*
*r*

*M*

*M* *E*

*S* *dr* *M*

*r* *r*

− + π − γ

### ∫

_{v}^{(5) }

Taken near the threshold γ ≈0, the transparency e^{−}^{2 /}^{S}^{}
gives the width Γ* _{n}*|

*E*

*| exp ( 2− π α*

_{n}*Z*), where α= e /(

^{2}κ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 |*

_{F}*Zα*| > 1/ 2 there are oscillations of the Ohmic conductivity which have a characteristic form of Fano resonances centered at

*E*

*[12]. In this regime the conductivity exhibits peaks at the densities for which the Fermi energy*

_{n}*E*

*equals*

_{F}*E*

*. The peak position is highly sensitive to the potential strength*

_{n}*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

( ) = *Ze* ,( > ), ( ) = *Ze*0 ,( < )

*V 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

= =

*z* *z* *z* 2 *z*

*J* *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 r*

*a* *j* *b*

*r*

*b E* *V r*

*b* *j* *a*

*r*

+ ∆ −

′ − + +

− ∆ −

′ + − −

*v*
*v*

(8)

It is convenient to define the quantities ε= /*E* v* _{F}*,

= / _{F}

*m* ∆ v , and α α= * _{g}* / = e /κ

^{2}v

*κ.*

_{F}The discrete spectrum of Eqs. (8) exists for | | <ε *m*_{. In }
this case it is convenient to define

2 2

= , = 2 , = ( ),

2

= ( )

2

*u* *m* *ur a* *m* *g f*

*b* *m* *g f*

− ε ρ + ε −

− ε +

(9)

and rewrite Eqs. (8) in the region r > r_{0} as follows:

1 = 0,

2 2

1 = 0.

2 2

*g g* *Z* *f j Z* *m*

*u* *u*

*f* *f* *Z* *g j Z* *m*

*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* *Z*

*d g* *u* *g*

*d*

+ αε − + α

+ − + +

ρ ρ ρ

(11)

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 *W*_{k}_{,}ν( ),*z M*_{k}_{,}ν( )*z* at infinity,

, ( ) e * ^{ur}*(2 ) ,

*W*

_{µ ν}ρ

^{−}

*ur*

^{µ}

, ( ) (1 ) e (2 ) ,

1 2

*M*_{µ ν} Γ + ν *ur* *ur* ^{−µ}
ρ Γ − µ + ν

*r*→ ∞, (13)

we find that the regularity condition at infinity requires 2= 0

*C* . Then the first equation in (10) gives the following
solution for the *f* component in the region II (*r r*> _{0}):

1 1 ,

2

= ( ).

*II* *Z*

*u*

*f* *C j Z* *m* *W*

*u* ^{− + α ν}^{ε}

− α ρ

(14)

Solutions in the region I (*r r*< _{0}) are easily obtained

2 2

1 | 1/2|

= 0 ,

*I* *j* *Z*

*b* *A rJ* *r* *m*

+ *r*

α

ε + −

(15)

2 2

1 0 | 1/2|

0 0

= sgn ( ) / ,

*I* *Z r m*/ *j* *Z*

*a* *A* *j* *rJ* *r* *m*

*Z r m* ^{−} *r*

ε + α + ε + α −

ε + α −

(16)
where *A*_{1} is a constant and we took into account the infra-
red boundary condition which selects only regular solution
for *b** _{I}* and

*a*

*. Energy levels are determined through the continuity condition of the wave function at*

_{I}*r r*=

_{0},

=0 =0

= ,

*I* *II*

*I* *r r* *II* *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* *k*

*j* *W*

*u*

+ αεν

− + αεν

ρ +

α −

− ρ

| 1/2|

0

0 | 1/2|

/ ( )

= sgn ( ) ,

/ ( )

*j*
*j*

*Z* *r m* *J*

*k* *j* *m*

*u* *Z* *r m J*

+

−

ε + α − ρ + ε

ε + α + ρ (18)

= (*Z* *r*_{0})^{2} *m r*^{2 2}_{0}.

ρ α + ε −

We analyze this equation in the limit *r*_{0}→0_{ where we }
can use the asymptotical behavior of the Whittaker func-
tion at ρ →0,

1 1

2 2

, ( ) (1(2 ) ) (1( 2 ) ) .

2 2

*W*_{µ ν} Γ ν ^{−ν} Γ − ν ^{+ν}

ρ ρ + ρ

Γ − µ + ν Γ − µ − ν

(19)

In the limit *r*_{0}→0 Eq. (18) reduces to the following one,

0 2

( 2 ) 1 (2 ) =

(2 ) 1

*Z* *u* *ur*
*Z* *u*

ν

ε

Γ + ν − α

Γ − ν

ε Γ ν Γ − ν − α

0

0 0

( ) ( )

= ( ),

( ) ( )

*Z m* *Z m*

*j* *k j*

*u* *u* *O r*

*Z m* *Z m*

*j* *k j*

*u* *u*

α + ε α − ε

+ ν − + − ν −

− +

α + ε α − ε

− ν − + + ν −

(20)

where

| 1/2|

0 | 1/2|

( )

= sgn ( ) ( , ).

( )

*j*
*j*

*J* *Z*

*m* *m*

*k* *j* *Z j*

*u J* *Z* *u*

+

−

+ ε α + ε

≡ σ α

α ^{ (21) }

Equation (20) can be rewritten in more convenient form

0 2

( 2 ) 1 (2 ) =

(2 ) 1

*Z* *u* *ur*
*Z* *u*

ν

ε

Γ + ν − α

Γ − ν

Γ ν Γ − ν − α ε

( )

( , )

= ( ) ( , ).

*j* *Z mu* *j* *Z* *Z j*

*Z m* *j* *Z* *Z j*

*j* *u*

α − ε

− ν − + ν − ασ α

− + ν − α − ε − ν − ασ α ^{(22) }

In the limit *r*_{0}→0 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]),

2 2 1/2

, 2

= 0,1,2,3,..., > 0,

= 1 ,

= 1,2,3,..., < 0.

( )

*n j* *m* *Z* *n* *j*

*n* *j*

*n*

α −

ε +

ν +

(23)

The bound states for *n*≥1 are doubly degenerate,

, = ,

*n j* *n j*−

ε ε . The lowest energy level is given by

0, =1/2* _{j}* =

*m*1 (2

*Z*) .2

ε − α (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 ε_{n}_{,1/2} be-
come complex for *Z*α> 1/ 2. The unphysical complex
energies indicate that the Hamiltonian of the system is not
a self-adjoint operator for supercritical values *Z*α> 1/ 2
and should be extended to become a self-adjoint operator.

According to [2,3], nonzero *r*_{0} resolves this problem. For
*Z*α> *j*, ν is imaginary for certain *j* and for such *j* we

denote ν= ,*i*β β= *Z*^{2 2}α −*j*^{2}. For finite *r*_{0} discrete levels
also exist for *Z*α> 1/ 2. Their energy decreases with in-
creasing of *Zα* until they reach the lower continuum. The
behavior of lowest energy levels with *j*= 1/ 2 as functions
of the coupling Zα is shown in Fig. 1(a).

The critical charge *Z** _{c}* that corresponds to diving into
the continuum is obtained from Eq. (22) setting =ε −

*m*there and using the corollary of the Stirling formula:

( ) e2 ln ,

( ) ^{iy x}

*x 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,

ln (2*Z mr*0) = arg (*Z* (*Z j*, ) *j i* )

−β α ασ α − + β +

arg (1 2 )*i* *n n*, 0,1, ...

+ Γ − β + π = _{(26) }

It is not difficult to check that for *j*= 1/ 2 and = 1*n* the
critical coupling *Z** _{c}*α approaches the value 1/2 for

*mr*

_{0}→0. The dependence of the critical coupling

*Z*

*α*

_{c}_{ on }

*mr*

_{0}for

= 1/ 2

*j* is shown in Fig. 1(b).

*2.2.2. Gapped graphene, supercritical regime *
Let us analyze Eq. (18) in the supercritical case

> 1/ 2

*Z*α and show that there are resonant states for

| | >ε *m* (we define the gap > 0∆ ). The Whittaker function
, ( )

*W*_{µ ν} ρ _{ with }µ= 1/ 2+ αε*Z* /*u*_{, }ν_{=} _{j}^{2}−_{Z}^{2 2}α describes
bound states for | |<ε *m* which are situated on the first
physical sheet of the variable u and for which Re > 0*u*
(see, Eq. (14)). The quasistationary states are described by
the same function *W*_{µ ν}_{,} ( )ρ and are on the second unphysi-
cal sheet with Re < 0*u* . 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 *Z*^{2 2}α > *j*^{2} resonance states are deter-
mined by Eq. (18) for bound states where ν is replaced by

=*i*

ν β. We will consider the states with *j*= 1/ 2 which
correspond to the *nS*_{1/2}-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* *W* *k*

*u*

*Z* *r m J*
*k* *m*

*u* *Z* *r m J*

*Z* *r* *m r*

+ αε β

− + αε β

ρ +

α −

− ρ

ε + α − ρ

+ ε

ε + α + ρ

ρ α + ε −

(28)

The analytical results can be obtained for the near-
critical values of Z when *Z*α −1/ 2 1 . We assume that

| 2*ur*0| 1 , then using the asymptotic of the Whittaker
function, we find

0 2

(1 2 ) 1

(2 ) =

(1 2 ) 1

*i*

*i* *Z*

*i* *u*

*ur* *i* *i* *Z*

*u*

β

αε

Γ + β −

Γ − β

αε Γ + β Γ − β −

1 0 1 0

( )

1 ( ) 1

2 ( )

= 12 ( ) 1 ( ).

2 2 ( )

*Z m* *i* *Z* *J Z*

*i* *u* *J Z*

*Z m* *J Z*

*i* *i* *Z*

*u* *J Z*

α − ε + β − α α

− β − α

α − ε α

+ β − − β − α

α

(29)

Expanding Eq. (29) in the near critical region in powers of

= *Z*2 2 1/ 4

β α − , 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 *Z*α< 1/ 2); black solid lines are numerical solution
for *j*= 1/ 2, *mr*_{0}= 0.01; black dashed line is a numerical solu-
tions for = 1/ 2,*j* − *mr*_{0}= 0.01 (a). The critical coupling as
a function of *mr*_{0} for the 1S_{1/2} level (b).

2 2 2 0

0 0 1

(1/ 2)

( 2 ) = 1 4

(1/ 2) (1/ 2)

*i* *J*

*i* *m r* *i*

*J* *J*

β

− ε − + + β − +

2 2

1 1

(1) 1

2 2

1
*i*

*m* *i* *m*

*m*

ε

+ Ψ − Ψ − ε − − + ε −ε + . (30)

Here Ψ( )*x* is the psi-function and we put *u*=− ε −*i* ^{2} *m*^{2}
where Im ε −^{2} *m*^{2} < 0 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 *Z** _{c}*α

_{ for }

*j*= 1/ 2 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)
*Z**c*

*cmr*
*c* *J*

*J* *J*

α + π

− Ψ − ≈

−

(31)

Note that the dependence of the critical coupling on *mr*_{0} is
quite similar to that in the strongly coupled QED [34,48].

For *Z Z*> * _{c}*, using Eq. (30), we find the following reso-
nant states:

3 / 2 3

= 1 e , = ,

8 * ^{b}* 8

_{c}

^{c}*m* *b i* π ^{−π} *b* πβ − β

ε − + + ββ (32)

where β* _{c}* = (

*Z*

*α −)*

_{c}^{2}1/ 4. Like in QED [49] the imagi- nary part of energy of these resonant states vanishes expo- nentially as

*Z*→

*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*

_{c}*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

*V**'*

*a r* *dr* *r*

*r* *m V*

− χ

ε + −

### ∫

^{}

_{}

^{(33) }

Here

2 2

2( ) = 2( ( )), = , = ,

2 _{F}

*m* *V*

*k r* −*U r* ε − *V*

v (34)

and we represent the effective potential as the sum of two
terms *U U U*= _{1}+ _{2}_{, where }*U*_{1} is the effective potential for
the Klein–Gordon equation and *U*_{2} takes into account the
spin dependent effects,

2

1 2

( 1)

= ,

2 2

*V* *j j*

*U* *V*

*r*

ε − + − (35)

2

2=1 3 2 .

4 2 ( )

*''* *'* *'*

*V* *V* *jV*

*U* *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 (*Z*→*Z** _{c}*,

*j*= 1/ 2, 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, = 0*m* .
Writing = | | eε ε ^{i}^{γ} Eq. (30) takes the form

ln (2 | | )0

*r* *i* 2π
ε + γ −

0

0 1

(1/ 2) 1 1

2 (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

*n* *ar* *i* *n*

*Z*

− γ π

ε −

α −

01 2 2

= (1.18 0.17 ) exp , = 1, 2, ,

1/ 4

*i r* *n* *n*

*Z*

− π

− + −

α −

(38) where

= 1 coth 3.28,

2 2

π π

γ + ≈

0

0 1

2 (1/ 2)

= exp1 2 (1) 1 Re 1 1.19.

2 (1/ 2) (1/ 2) 2

*J* *i*

*a* *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” *e*^{2}/*v** _{F}* ≈2.2 in gra-
phene, an instability could potentially appear already for
the charge = 1

*Z*. However, in the analysis above we did

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 *Z*_{eff} is such
that the impurity with bare charge = 1*Z* remains subcriti-
cal, *Z*_{eff}e /(^{2} κv* _{F}*) < 1/ 2, for any coupling

*e*

^{2}/(κv

*), while impurities with higher*

_{F}*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 *m*^{2 2}/ε ,

2 2

(0) 2

= 1 (0.29 0.23 ) ,

2 ^{n}

*m* *m m* *i*

ε − ε ε − ε + ε − *n*= 1, 2,. (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 = 0

*m*. 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 }

*m*

*and easily find the first two terms*

^{k}(0) 2

= (0) (0.24 0.20 ).

| |

*n* *n*

*n*

*m* *m* *i*

ε ε − + +

ε (41)

Since Imε^{(0)}* _{n}* < 0, 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

*Z*

*α(*

_{c}*mr*

_{0}), 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.62 in the static random-phase approximation [21] and α*

_{c}*= 0.92*

_{c}_{ 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

*Z*

*α in the problem of the Cou-*

_{c}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 NH_{3}, 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 *T* < 10 K. 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 (*Z Z*> * _{c}*).

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 *Z** _{c}*=

*v*

*/ 2e*

_{F}^{2}0.25. By tuning the graphene Fermi level

*E*

*via electrostatic gating the atomic collapse behavior was observed.*

_{F}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

*( , ) is given by Eq. (64) below. The curves in Fig. 2 show the differential conductance*

**D E r***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

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, *Z Z*/ * _{c}* 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 *Z Z*/ * _{c}* 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

*Z Z*/

*= 0.6 0.3± . This is in agreement with the value*

_{c}*Z Z*/

*= 0.5 0.1± 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].*

_{c}**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, *Z** _{c}*α 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 (*V** _{g}*= 30 V− ) 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

*Z 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.*

_{c}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 = 0*n* , *j*= 1/ 2− comes close to the
highest energy state of the level = 1*n* − . 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 *n*≤ −1.

*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 **p**→**p**=− ∇ +*i* ^{e}_{c}**A**, where

< 0
*e*

− 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 *r*_{0} 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

( ) = *Z*e

*V* − *r* *r*

κ +

**r** . (42)

It is convenient to introduce the magnetic length

= / | |

*l**B* *c eB* and the dimensionless quantity ζ= e /(*Z* ^{2} κv* _{F}*)
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

*B* *F*

*B* *F*

*j* *r* *E* *V r*

*a* *a* *a* *b*

*r* *l*

*j* *r* *E* *V r*

*b* *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 ( )

[*E* *V r*( )] ( ) =*r* *b r* .

− ξ∆ − χ *r* _{(44) }

We obtain the following Schrödinger-like equation:

( )*r U r r*( ) ( ) = ( ),*r*

−χ′′ + χ χ (45)

where

2 2

( * _{F}*)2

*E* − ∆

= *v*

(46)

and the effective potential, *U U U*= _{1}+ _{2}_{, reads }

2

1 2 2 4 2

(2 ) ( 1) 1/ 2

= ,

( * _{F}*) 4

_{B}

_{B}*V E V* *j j* *r* *j*

*U* *r* *l* *l*

− + −

+ + +

v ^{ (47) }

2

2=1 3 2 2 .

2 2 2* _{B}* ( )

*V* *V* *j* *r* *V*

*U* *E* *V* *E* *V* *r* *l* *E* *V*

′′ ′ ′

− ξ∆ − + − ξ∆ − − + − ξ∆ −

(48)
We plot the effective potential *U r*( ) near the *K*_{−} point
( = 1ξ − _{) for =}* _{E}* −∆ and

*j*= 1/ 2− 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 = 0*B* .

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 *Z*α →0, we can find it using per-
turbation theory. For *Z*α= 0, 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 *E** _{nj}*.
The energy downshift is largest for

*E*

_{0}

*and diminishes with increasing –j. For the = 0*

_{j}*n*level with

*E*

_{0}

^{(0)}

*=*

_{j}*E*

_{0}

^{(0)}=∆ 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/ 2*j* − .

2/42 2 /2

2

( 1) 0

( , ) = e ,

2 ! e

2

*M* *r* *lB* *M*

*M* *iM*

*B* *B*

*r* *l* *M* *r*

*l*

−

− φ

−

Ψ φ

π (49)

where *M* = (− +*j* 1/ 2) = 0, 1, 2, ... is the orbital quantum
number. Energy corrections ε* _{M}* of perturbed states of the
Landau level

*E*

_{0}

^{(0)}=∆ are found from the secular equation

|ε −*V** _{M M}*1 2 | = 0, (50)
where

*V*

*M M*

_{1 2}is a matrix element of the potential on states (49). Since

*V*

_{M M}_{1 2}is a diagonal matrix, we easily obtain

2 2 1 2/2

2 2

0 0

= = =

!2

*M*

*M* *MM* *M*

*B*

*Ze* *d* *e*

*V* *M* *l*

∞ ρρ + −ρ

ε −

κ

### ∫

ρ + ρ### ( )

2 2 1 2

0 0

= 1 1 ,3 / 2 ; / 2 ,

2

*M* *M*
*B*

*Ze* *M* *M*

*l*

+

− + ρ Ψ + + ρ

κ (51)

where Ψ( , ; )*a c z* is the confluent hypergeometric function,
0= /*r l*0 _{B}

ρ . For small *Zα*1 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*=∆,

(0) 2

0 = 0 .

*M* *M* 2

*B*

*E* *E* *Ze*

*l* *M*

+ ε ∆ −

κ (53)

The largest correction by modulus ε_{0}=− α*Z* v* _{F}* π/

*l*

*2 is for the state with*

_{B}*M*= 0. 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 *E E*= _{0}^{(0)}+ ε_{0}=−∆ 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, *Z*α1 and ∆*l** _{B}* 1, this gives

= 2 2 * ^{B}*.

*c* *F*

*Z* ∆*l*

α πv (54)

Clearly, this critical charge tends to zero as ∆ →0, while
the state with *M* = 0 of the zero Landau level moves be-
low zero energy for any small impurity charge (its energy
is ε_{0}=− α*Z* v* _{F}* π/

*l*

*2). 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).*

_{B}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 *r*→0 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 = 1

*m j*− − −*M* . Then, for *m*≥0(*j*≥1/ 2) the upper
component of (7) dominates and the leading behavior is

( ) * ^{m}* 1

*a r* *r* ^{+}, for < 0*m* (*j*< 1/ 2) the lower component
dominates with ( )*b r* *r*^{−}* ^{m}* (see Ref. 67).

The numerical integration of Eq. (43) proceeds as fol-
lows. We take a “shot” from = 0*r* 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 = 10*B* 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*Ψ, = 0*m* and = 1*n* − , = 0*m* levels, dashed lines correspond to
anticrossing of the = 1*n* , = 1*m* − and = 1*n* − , = 1*m* − levels. In both panels the regularization parameter is chosen as *r*_{0}= 0.5 nm (b).