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

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

(Review Article)

E.V. Gorbar

, V.P. Gusynin

, and O.O. Sobol

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

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

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 )m_e ²/ (2m r_e ²) 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ⁿ with n≥2. In fact, the Schrö- dinger equation with the potential V r( ) =−β/r² provides the canonical textbook example of the fall-to-center in quan- tum mechanics [1], which takes place for β>²/ (8 )m_e when the energy spectrum is not bounded from below. If the interaction potential is regularized at some distance r₀, then the electron wave function of the ground state is local- ized in the region of the radius r₀ which shrinks to the origin as r₀→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 /(² v_F) 2.2≈ , where 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 /² 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_cα 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 z} ( ),

H p ξ v σp+ ξ∆σ +V r ₍₁₎ where v_F is the Fermi velocity of graphene, p=−i∇ 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 Θ=is_{2 1}σK_:

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

H ⁻ H

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

where s₂ 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 ² c.

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

= _F | | ( )

H v p +V r , where V r( ) =−Ze / ( )² κr and κ is a di- electric constant. The collapsing trajectories with angular momenta M M< _c = e / (Z ² κv_F) are separated from non- falling trajectories by a centrifugal barrier. This is mani- fested in the expression for the radial momentum square

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₁< <r r₂, r_1,2= ( e /Z ² κ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₀ is a regularization parameter, which is of order of lattice spacing. Evaluating the integral with logarithmic accuracy, we obtain γln ( e / (Z ² κr E₀| |)) =πn, where

2 2 1/2

(M_c M )

γ ≡ − , which gives the quasi-Rydberg states

2 /

e e0 ⁿ , > 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 Imp_r and obtain the tunnel- ing action

( )

2 2 2

2 1

= = .

r c

F c r

M

M E

S dr M

r r

 

− +  π − γ

 

∫

Taken near the threshold γ ≈0, the transparency e^{−2 /S } gives the width Γ_n|E_n| exp ( 2− π αZ ), where α= e /(² κv_F) 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 |Zα| > 1/ 2 there are oscillations of the Ohmic conductivity which have a characteristic form of Fano resonances centered at E_n [12]. In this regime the conductivity exhibits peaks at the densities for which the Fermi energy E_F equals E_n. 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

( ) = Ze ,( > ), ( ) = Ze0 ,( < )

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

φ − φ +

 

 

Ψ  

  ⁽⁷⁾

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 /κ ² 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₀ 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> ₀):

1 1 ,

2

= ( ).

II Z

u

f C j Z m W

u ^{− + α νε}

 − α  ρ

 

  (14)

Solutions in the region I (r r< ₀) 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₁ is a constant and we took into account the infra- red boundary condition which selects only regular solution for b_I and a_I. Energy levels are determined through the continuity condition of the wave function at r r= ₀,

=0 =0

= ,

I II

I r r II r r

b b

a a ⁽¹⁷⁾

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₀)² m r^{2 2}₀.

ρ α + ε −

We analyze this equation in the limit r₀→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 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

+

−

+ ε α + ε

≡ σ α

α ⁽²¹⁾

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

α − ε

− ν − + ν − ασ α

− + ν − α − ε − ν − ασ α ⁽²²⁾

In the limit r₀→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 (2Z ) .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₀ resolves this problem. For Zα> j, ν is imaginary for certain j and for such j we

denote ν= ,iβ β= Z^{2 2}α −j². For finite r₀ 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

− β α β − + ασ α Γ − β

− β − + ασ α Γ + β ⁽²⁵⁾ or,

ln (2Z mr0) = arg (Z (Z j, ) j i )

−β α ασ α − + β +

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

+ Γ − β + π = ₍₂₆₎

It is not difficult to check that for j= 1/ 2 and = 1n the critical coupling Z_cα approaches the value 1/2 for mr₀→0. The dependence of the critical coupling Z_cα _on mr₀ 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}²−_Z^{2 2}α 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 Z^{2 2}α > j² 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

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

= Z2 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.01; black dashed line is a numerical solu- tions for = 1/ 2,j − mr₀= 0.01 (a). The critical coupling as a function of mr₀ 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 ² m² where Im ε −² m² < 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) Zc

cmr c J

J J

α + π

 

− Ψ − ≈

 − 

 



(31)

Note that the dependence of the critical coupling on mr₀ 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α −)² 1/ 4. Like in QED [49] the imagi- nary part of energy of these resonant states vanishes expo- nentially as Z→Z_c. 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

V'

a r dr r

r m V

   

− χ

  ε + −  

   



∫

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= ₁+ _{2, where} U₁ is the effective potential for the Klein–Gordon equation and U₂ 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

   

 +   + 

ε + − ε + −  ε + − 

 

  

   ⁽³⁶⁾

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, = 0m . 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²/v_F ≈2.2 in gra- phene, an instability could potentially appear already for the charge = 1Z . 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 = 1Z remains subcriti- cal, Z_effe /(² κv_F) < 1/ 2, for any coupling e²/(κv_F), 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 m^{2 2}/ε ,

2 2

(0) 2

= 1 (0.29 0.23 ) ,

2 ⁿ

m  m m i 

ε − ε ε  − ε + ε −  n= 1, 2,. (40) The resonant states with ε⁽⁰⁾_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

ε

∑

(0) 2

= (0) (0.24 0.20 ).

| |

n n

n

m m i

ε ε − + +

ε (41)

Since Imε⁽⁰⁾_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₀), 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 α_c = 1.62 in the static random-phase approximation [21] and α_c = 0.92 _{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_cα 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 NH₃, 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_F/ 2e²0.25. By tuning the graphene Fermi level E_F 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

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/ _c = 0.6 0.3± . This is in agreement with the value Z Z/ _c= 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].

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/ _c 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.

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 , j= 1/ 2− 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 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_cA, 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₀ 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

( ) = Ze

V − r r

κ +

r . (42)

It is convenient to introduce the magnetic length

= / | |

lB c eB and the dimensionless quantity ζ= e /(Z ² κ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 ₍₄₄₎

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= ₁+ _{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 ⁽⁴⁷⁾

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 = 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 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_0j and diminishes with increasing –j. For the = 0n level with E₀⁽⁰⁾_j =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 − .

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₀⁽⁰⁾=∆ are found from the secular equation

|ε −V_{M M}1 2 | = 0, (50) where VM M_{1 2}is a matrix element of the potential on states (49). Since V_{M M1 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 l0 _B

ρ . For small Zα1 we can use the unregularized Coulomb potential, then setting ρ₀ = 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 ε₀=− αZ v_F π/l_B 2 is for the state with 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= ₀⁽⁰⁾+ ε₀=−∆ 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 ε₀=− αZ v_F π/l_B 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).

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 < 0m (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 = 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 r₀= 0.5 nm (b).