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

We review the theoretical and experimental results connected with the electron states in two-dimensional Dirac systems paying a special attention to the atomic collapse in graphene. Two-electron bound states of a Coulomb impurity are considered too. A rather subtle role of a magnetic field in the su...

Full description

Saved in:

Bibliographic Details
Date:	2018
Main Authors:	Gorbar, E.V., Gusynin, V.P., Sobol, O.O.
Format:	Article
Language:	English
Published:	Фізико-технічний інститут низьких температур ім. Б.І. Вєркіна НАН України 2018
Series:	Физика низких температур
Subjects:	Обзор
Online Access:	https://nasplib.isofts.kiev.ua/handle/123456789/176115
Tags:	Add Tag No Tags, Be the first to tag this record!
Journal Title:	Digital Library of Periodicals of National Academy of Sciences of Ukraine
Cite this:	Electron states in the field of charged impurities in two-dimensional Dirac systems / E.V. Gorbar, V.P. Gusynin, O.O. Sobol // Физика низких температур. — 2018. — Т. 44, № 5. — С. 491-524. — Бібліогр.: 142 назв. — англ.

Institution

Digital Library of Periodicals of National Academy of Sciences of Ukraine

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

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

Institution

Similar Items