The critical transition of Coulomb impurities in gapped graphene

The effect of supercritical charge impurities in graphene is very similar to the supercritical atomic collapses in QED for Z>137, but with a much lower critical charge. In this sense graphene can be considered as a natural testing ground for the analysis of quantum field theory vacuum instabilities. We analyze the quantum transition from subcritical to supercritical charge regimes in gapped graphene in a common framework that preserves unitarity for any value of charge impurities. In the supercritical regime it is possible to introduce boundary conditions which control the singular behavior at the impurity. We show that for subcritical charges there are also non-trivial boundary conditions which are similar to those that appear in QED for nuclei in the intermediate regime 118<Z<137. We analyze the behavior of the energy levels associated to the different boundary conditions. In particular, we point out the existence of new bound states in the subcritical regime which include a negative energy bound state in the attractive Coulomb regime. A remarkable property is the continuity of the energy spectral flow under variation of the impurity charge even when jumping across the critical charge transition. We also remark that the energy levels of Hydrogenoid bound states at critical values of charge impurities act as focal points of the spectral flow.


Introduction
The stability of non-relativistic hydrogenoid atoms is one of the essential features that contributed to consolidate the quantum theory. However in relativistic quantum mechanics there is a critical value of the central point-like charge Ze from where on atom stability is lost [1]- [4]. This is one of the surprising consequences of relativistic invariance in QED. The phenomenon can be understood in a heuristic way as a falling to the center catastrophe. The critical value in QED is reached when the spectrum of bound states of Dirac equation becomes complex which occurs for Z > 137. In fact what happens in QED is that when one of the bound states reaches the negative continuum spectrum the vacuum becomes unstable, generating electron-positron pairs. The positron escapes to infinite and the electron screens the central charge. The phenomenon has attracted attention from a fundamental viewpoint because it suggests that could be experimentally tested by detecting an excess of positron in the collisions of heavy nuclei [5] [6].
The instability of the atom for supercritical charges has also inspired a new mechanism of quark confinement in QCD [7]- [11]. The running of the strong coupling under the renormalization group flow in QCD leads to large values of the effective charge of the quarks, which reaches very fast supercritical values in the infrared. The instability of the vacuum generates a transition from the perturbative Coulomb regime at short distances to a confinement regime in quarks interactions at large distances [12]- [15].
The discovery of graphene [16] opened a new window for the analysis of this phenomenon [17]. In that case a similar phenomenon occurs in the presence of charged impurities, but with a much lower critical charge. In graphene the instability yields to a screening of the charge impurity, and the phenomenon has been recently experimentally observed [18,19]. Motivated by this new physical effect we review the main features of this phenomenon and shed some light in some of its more paradoxical aspects. We address the problem from a viewpoint where quantum unitarity is never lost no matter the strength of charge impurities. In fact we show that the formal analytic continauation of the bounded energy levels of the Coulomb problem into complex values does not mean a loss of Hermiticity of the corresponding effective Hamiltonian. It only shows the existence of non-trivial spectral densities in the continuum spectrum which correspond to the existence of resonances in scattering processes [18,19].
In order to clarify this issue we analyze in graphene the transition from the subcritical regime to the supercritical one by increasing the values of impurity charges. The results show a continuous behavior of the corresponding energy levels, although the spectral flow is very peculiar: energy levels of Hydrogenoid spectrum in the critical regime are focal points of the spectra of subcritical and supercritical regimes. The peculiar behavior of the supercritical regime is reflected by the increasing number of energy levels inside the energy gap, but the continuity of the spectral flow is always preserved along the transitions between the different spectral regimes. Vacuum instability of the corresponding quantum field theory is pointed out by the crossing of the E = −m energy level of the Dirac sea continuum by some eingenvalues of the Dirac Hamiltonian which implies the appearance of pair particle-antiparticle creation mechanism that leads to the screening of the charge impurity.
The analysis of the problem is based in a novel method of dealing with selfadjoint extensions of the Dirac Hamiltonian. In that formalism all cases are approached in an unified and global way that allows to follow the spectral flow of the different (weak-strong) regimes in a smooth way. The analysis can be extended to any space dimension, e.g see [20] for the three dimensional case.
In the Section 2 we analyze the unitarity problem of the Dirac Hamiltonian in a Coulomb background. The problem is solved by using the theory of self-adjoint extensions which regularize the singularities associated to the Coulomb potential. The selfadjoint Hamiltonians are classified in different regimes according to the value of impurity charges. In Section 3 we calculate the bound states energy spectrum of the Coulomb Hamiltonian in the different regimes. A particular attention is paid to the special cases of Hydrogen and meta-Hydrogen spectra (see [21][22]), The spectral flow of the bound states spectrum is analysed in Section 4, where we also study the analytic properties of this flow in the different subcritical and critical regimes. Finally, the analysis of the results and conclusions is carried out in Section 5.

Charged impurities in graphene
Graphene is a two dimensional layer of carbon atoms arranged on a honeycomb lattice of hexagons. The magic of the hexagonal honeycomb structure of graphene leads to a spectral structure in the first Brillouin zone with two contact points K and K (Dirac points) between electronic bands. In a neigbourghood of any of these two points 1 the spectrum of unbounded electrons is well described in terms of a massless Dirac Hamiltonian [23] where σ i , i=1,2,3, are the Pauli matrices and v F is the velocity of the electrons at the Fermi surface, which for suspended graphene is about 300 times smaller than the speed of light in vacuum. This behavior also holds for graphene in a substrate of SiO 2 with a slight modification of v F . Although natural graphene behaves like a semi-metal with no spectral gap, for electronic applications it is convenient to open a gap between the bands to reach a semiconductor regime. This behavior can be attained by different methods, either by introducing some disorder or by epitaxially grow graphene on a SiC substrate [24]. In that case the effective Hamiltonian (10) becomes a massive Dirac Hamiltonian where m the effective mass of the gap.
In the presence of a charged Coulomb impurity the effective electronic Hamiltonian becomes where is the effective charge of the impurity, the effective dielectric constant of the graphene sheet, r = x 2 + y 2 and the electronic effective speed factor v F has been absorbed by rescaling of The values of α depend on the substrate where the graphene sheet is grown. For instance, α ≈ 2 for vacuum, α 1 for SiO 2 and α 0.35 for SiC. The presence of a charge impurity with strong Coulomb interactions generate remarkable effects in the spectroscopic and transport properties. The physics of the effective theory is quite similar to that of relativistic atomic physics where the presence of instabilities is rather well know [1]- [4]. In any case there is a renewed interest on the theoretical and experimental studies on the Coulomb potential supercritical instabilities [25]- [39]. The main difference with respect to the 3D analogue (Hydrogen-like atoms) is that the value of the supercritical charge is much smaller α = 1 2 137. Although the single particle approach to the Coulomb problem constitutes the first step in addressing nontrivial features of the full-fledged many-body interacting theory, most of the phenomenology of graphene physics can be explained from this simplified approach. 1 For simplicity we consider only the K Dirac point and = 1

2D Dirac Hamiltonian in a Coulomb background
The presence of a singularity at the origin of the Coulomb potential requires the use of some ultraviolet renormalization mechanism. For such a reason it is convenient to introduce an ultraviolet cut-off r < > 0 around that singular point r = 0 and later on take the appropriate limit to extend the Hamiltonian to the whole space R 2 \{0} [1, 40,41]. The only physical requirement is unitarity of time evolution, which is equivalent to the self-adjointness of the Hamiltonian defined in such limit. If we exclude from the physical space a disk D(r < ) = {x ∈ R 2 ; x < r < } of radius r < around the origin, the most general boundary conditions can be given by [42,43,44] (1 + n /)ψ( in terms of a unitary operator U (r < ) defined on the boundary values of spinors ψ(r < ) ∈ L 2 (S 1 r< , C 4 ), wheren denotes the normal vector to the circumference S 1 r< = {x ∈ R 2 ; x < r < }. Using polar coordinates r and θ a general spinor ψ can be expanded as in terms of orthogonal eigenfunctions of the total angular momentum J z = L z + S z = −i ∂ ∂φ + 1 2 σ 3 , with semi-integer eigenvalues j = l + 1/2. The space of spinors can be then decomposed as orthogonal sum of subspaces with fixed total angular momentum j ∈ Z + 1 2 : where ψ j is a spinor of the form which belongs to the subspace of total angular momentum j, In order to preserve the SO(2) rotation symmetry in the regularized theory, the unitary operator U (r < ) fixing the boundary condition has to be diagonal in the angular momentum decomposition (5), i.e. on each subspace of fixed angular momentum j the unitary operator U (r < ) reduces to a single phase e 2i β j . Thus, the boundary condition (4) becomes More explicitly, where F (r < ) and G(r < ) are real functions. The removal of the UV regularization requires to take the limit r < → 0 which implies the choice of an appropriate series of boundary conditions U (r < ). The optimal choice of boundary conditions U (r < ) that guarantees the convergence of the UV limit is given by the flow driven by asymptotic zero modes. Near the impurites asymptotic zero modes are solutions of the equation They will play a fundamental role in the renomalization of the singularity introduced by the impurities as they do in the three-dimensional case of Hydrogenoid atoms [20][45] [46]. The key observation is that in the vicinity of the inpurity 0 < r r < any solutions of the Coulomb-Dirac equation Hψ = Eψ behaves as an asymptotic zero mode. Thus, all the spinors in the domain of the Hamiltonian must behave near the singularity as zero modes of (10).
For any choice of boundary condition β j 0 at a given cut-off r 0 1 there is a unique asymptotic zero mode (F j 0 , G j 0 ) satisfying the equation If the two components of the asymptotic zero mode (F j 0 , G j 0 ) are L 2 normalizable in a neigbourghood of the singularity, i.e. F j 0 , G j 0 ∈ L 2 (D(r < ), C), then the flow of boundary conditions β j r< (r < ∈ (0, r 0 )) given by . (12) defines in the limit r < → 0 a selfadjoint extension of the Dirac Hamiltonian (3). The domain of the Hamiltonian is expanded by the spinors (F j 0 , G j 0 ) which satisfy In other terms, once the cut-off r 0 is fixed, we can associate to each boundary condition parametrized by β j 0 a unique asymptotic zero mode satisfying (11). On the other way, given an asymptotic zero mode, the relation (12) defines for each r < ∈ (0, r 0 ) a boundary phase β j r< in a unique way.
For some values of the impurity charge not all the boundary conditions β j 0 give rise to normalizable asymptotic zero modes. Such boundary conditions do not lead by the procedure described above to a well defined selfadjoint Dirac Hamiltonian. However, as we shall see later on, it is always possible to find an alternative boundary condition β j 0 for the same value of impurity charge whose zero mode is normalizable and leads to a well defined selfadjoint Dirac Hamiltonian 2 .
By this method we have replaced the convergent flow of UV cut-off boundary conditions just by the choice of a simple asymptotic boundary condition (13). The boundary condition flow is then defined in this way: the initial cut-off phase β j 0 defines an asymptotic zero mode (F j 0 , G j 0 ), and the boundary phases β j (r < ) run with the cut-off while keeping fixed the zero mode, converging to a well defined boundary condition when the cut-off is removed.
In summary, the boundary condition of the Dirac Hamiltonian in a Coulomb background is defined by the choice of one of these two equivalent boundary data: either a unitary matrix U (r 0 ) of the form (7) acting on the functions of the boundary of a small cut-off disk of radius r 0 1 or a normalizable asymptotic zero mode (F j 0 , G j 0 ). The connection between the two choices is given by equation (4). Moreover, any boundary condition that leads to selfadjoint extension of the Hamiltonian (3) is obtained by this method.

Boundary conditions for different regimes
The subspace of asymptotic zero modes (F j 0 , G j 0 ) satisfying the boundary condition (9) for a given angular momentum j ∈ Z + 1 2 depends on the value of the charge α of the impurity, in a similar way as in the three-dimensional analogue case [20,45,46] To find the asymptotic zero modes of the Hamiltonian (3) we have to look only at leading terms asymptotic expansion around the impurity. Using the expansion (6) is easy to show that they satisfy the following coupled equations Searching for solutions of the form F j 0 (r) = r s and G j 0 (r) = C r s we find two independent solutions where ν = j 2 − α 2 . For α 2 = j 2 the two solutions degenerate, but in this case the logarithmic corrections give rise also to two independent solutions of the form where Λ is a parameter with dimensions of momentum [L] −1 .
Notice that the value j 2 = α 2 is critical: when α 2 < j 2 the parameter ν is real, while for α 2 > j 2 it is purely imaginary. Thus, depending on the strength of the charge impurity there are three different regimes where to impose the boundary conditions. a) Regular regime: This regime is never reached in the lowest angular momentum states j = ±1/2.
In this case ν is a real parameter and one of the two asymptotic zero modes solutions is not normalizable in a neigbourghood of the origen D(r 0 ). Indeed, the asymptotic zero mode is not square integrable in D(r 0 ). Thus, we are left with only one asymptotic behaviour given by the normalizable zero mode which strongly constrains the boundary condition (9), In particular, the parameter β j 0 is completely fixed, independently from r 0 , by This means that there is a unique self adjoint extension of the Hamiltonian (3). The boundary condition (13) becomes: lim In this regime both solutions are normalizable, thus the most general asymptotic zero mode is a linear combination of the two solutions (19) (20). The choice of β j 0 ∈ [0, π) fixes that linear combination in a unique way, up to a global constant.
where the parameter θ ∈ [0, π) of the linear combination is given by Thus, the boundary condition (13) becomes: In this case the most general asymptotic zero mode is where the parameter θ ∈ [0, π) can be related to the phase β j 0 ∈ [0, π) of the boundary condition (9) imposed at S 1 .
The corresponding boundary condition (13) of the Dirac operator becomes: In this case the value of ν becomes imaginary and both asymptotic zero modes are normalizable. A general zero mode solution is of the form where the parameter θ ∈ [0, π) is fixed by the phase β j 0 ∈ [0, π) of the boundary condition (9) imposed at S 1 The asymptotic boundary condition (13) in this case reads Notice that in any of the above regimes the Hamiltonian (3) is a selfadjoint operator. From now on we will parametrize the boundary conditions by θ ∈ (0, π) keeping in mind its relations with β j 0 and r 0 .

Bound states and energy levels
Once we have shown that the Dirac Hamiltonian (3) is a selfadjoint operator it is possible to analyze its spectrum by finding the energy levels The eigenvalue problem can be reduced, by using the ansatz (6) for each subspace of fixed angular momentum j, to solve the pair of coupled differential equations Let us now introduce two radial functions a(r) and b(r) defined by and use the notation = √ m 2 − E 2 , x = 2 r. The coupled equations satisfied by a(x) and b(x) become which can be easily decoupled. Indeed it is obvious to realize that and then The general solution of (35) can be expressed in terms of Whittaker functionsW and M [47] a(x) where A and B are constants. In the same way we have that Thus, the general solution of (29) and (30) is given by The asymptotic behavior of these solutions is strongly dependent on the regime of charge impurities.
If α 2 = j 2 the asymptotic behavior can be derived from the behavior of the Whittaker functions W and M for r 1, whereas for α 2 = j 2 : where ψ(x) = Γ (x)/Γ(x) is the digamma function and γ the Euler's constant. The spectrum of energy levels is also strongly dependent on the regime of charges. Let us focus on the discrete energy spectrum which correspond to bound states.

Regular regime
The boundary conditions (21) can only be satisfied if the constant A of the general solution (39) vanishes. On the other hand bound state spinors ψ have to be L 2 (R 3 , C 2 )-normalizable which implies that it must to decay at infinity. Thus, the asymptotic behaviour at r 1 of (29) and (30) must be of the form F j (r) ∼ = e − r , G j (r) ∼ = e − r . The implies that the spinors should look like where f (r) and g(r) are two radial functions that are polynomially bounded at infinity. This requirement is satisfied when the expressions This is the well known Hydrogenoid atom spectrum of bound states. For α 2 > j 2 − 1 4 the boundary conditions (23) and (25) for θ = 0 and θ = π 2 (we will analyze these two exceptional cases later separately), and (27) (for any value of θ) are satisfied only if the parameter B of the general solution (39) vanishes B = 0. In that case only terms involving the Whittaker function W survive, which implies that they automatically decays exponentially at infinity.
In this sense the exponential decay e − r of bound states means that they are localized around the impurity charge and thus behave as edge states in topological insulators [42,43].
Notice that in the massless limit the exponential decay e − r becomes a pure phase factor e −iEr and the corresponding solution is not localized and in fact belongs to the continuum energy spectrum.

Subcritical regime j
In that case the boundary conditions (23) are satisfied only if The solution of (45) gives the spectrum E II n (θ) of bound states in the subcritical regime.
3.3 Critical regime α 2 = j 2 , for θ = 0 and θ = π In this case the spectral condition derived from the boundary conditions (25) is where E = E/ , m = m/ and Λ = Λ/ . The solutions of equation (46) give an infinite sequence E III n (θ) of discrete energy levels.

Supercritical regime α 2 > j 2
In this case the infinite set of energy levels E IV n (θ) n ∈ Z is given by the spectral condition In the subcritical regime j 2 − 1 4 < α 2 < j 2 (3.2), we have two special cases: θ = 0 and θ = π 2 where the asymptotic zero modes that defines the boundary conditions are reduced to one of the two different asymptotic behaviors near the origin.
3.5 Subcritical regime j 2 − 1 4 < α 2 < j 2 with θ = 0 (Hydrogenoid atom) The boundary conditions reduce in this case those of the regular regime (3.1) , and then the spectrum is the same as in (44), i.e the bound states spectrum is the same as the Hydrogenoid atom E H n .
These boundary conditions can be satisfied by setting A = 0 in the general solutions and making the replacement ν → −ν. Using the same techniques as in the regular case, we get the analytic spectrum of bound states and with n = 1, 2, 3, .... The above bound states are known meta-Hydrogenoid states 3 . The meta-Hydrogenoid bounded spectrum is very similar to the Hydrogenoid spectrum. The only difference is a sign in the second radical. As in the Hydrogen case for α 2 > j 2 the spectral formulae (49) (50) becomes complex, but as we have already remarked the real spectrum is given by (47).

Critical regime
In the critical regime for boundary conditions with θ = 0 the Hydrogenoid and meta-Hydrogenoid spectra do coincide. They are defined by (44) with the only difference that E 0 = 0 for α 2 = j 2 and j > 0.
The Hydrogenoid and meta-Hydrogenoid spectra are not defined for α 2 > j 2 .

Spectral flows of bound states
The problem which inspired the Gribov approach to confinement is the fact that the energies of the bound states given by E H n become complex for α 2 > j 2 . To better understand that mechanism let us analyse the flow of the bound state spectrum by continuoulsy increasing the charge α of the impurity or varying the boundary conditions.

Spectral flow and boundary conditions
It is interesting to analyze the flow of the spectrum as we change the boundary condition parameter θ. The continuous flow of the spectrum defined by the change of the parameter θ characterizing the boundary conditions in the subcritical regime j 2 − 1/4 < α 2 < j 2 is displayed in Figure 1. There we plot for j = 3 2 the θ dependence of the lowest energy bound states in this regime.
Notice that in the limits θ = 0, θ = π/2, θ = π we recover the Hydrogenoid and meta-Hydrogenoid spectrum, i.e. The continuity of the flow should be obvious from the fact that the boundary conditions (23) reduce to the boundary conditions of the Hydrogenoid and meta-Hydrogenoid spectra in these limits. What is more surprising is the fact that the spectral flow is not periodic, i.e. there is an spectral asymmetry. The spectrum is periodic, i.e. it is the same at θ and θ + π, but the flow shifts the energy levels by one unit in each cycle from θ = 0 to θ = π. In general, for fixed angular momentum and charge we have for any integer k ∈ Z. The behaviour of the spectral flow recalls the pumping mechanism exhibited by edge states in topological insulators [48,49].
Another interesting property of the spectral flow is its monotoncity, i.e. E II (θ) < E II (θ ) if θ < θ . In particular this implies the standard sandwich inequalities between the Hydrogenoid and meta-Hydrogenoid energy levels Notice that there is a bound state emerging from the continuum E ≤ −m at a value of θ close to θ = 0. The behaviour of the flow is the same for positive j > 0 and negative angular momentum j < 0, except for the absence of zero levels (n=0) for the Hydrogenoid and meta-Hydrogenoid energy levels for j < 0 . Thus the sandwich inequalities in the negative case are Another interesting feature of the subcritical regime is that from (44) and (50) it follows that for n > 0 the spectra E H n and E h n with j > 0 and j < 0 are degenerate. The boundary condition (23) for θ = kπ and θ = 2k+1 2 breaks this degeneracy and creates a gap between the energies corresponding to j > 0 and j < 0. The situation is described in Figure 2 for Zα = 1.45. For θ = 0 we have the energy corresponding to n = 1 of the Hydrogen spectrum E H 1 which is degenerate for j = ± 3 2 . As we increase the parameter θ a gap appears between the states j = ± 3 2 . The energy of the state j = − 3 2 becomes lower than the one corresponding to j = 3 2 . The gap disappears again for θ = π 2 , where we have the again a degenerate energy level corresponding to n = 2 of the meta-Hydrogenoid spectrum E h 2 . If we increase the boundary condition parameter θ the gap reappears again. This time with the energy corresponding to j = 3 2 lower than the one corresponding to j = − 3 2 . Finally, for θ = π the two energy levels become again degenerate at the level n = 2 of E H 2 .
In the critical regime, α 2 = j 2 , as we have anticipated, the Hydrogenoid E III n (0) and meta-Hydrogenoid E III n (π) spectra do coincide and are given by (44) with the only difference that for α 2 = j 2 and j > 0, E 0 = 0. Once more this fact we can be understood in a simpler way, just by looking at the corresponding boundary conditions. Analyzing how the spectrum E III n (θ) changes with θ, we find that the correspondence in this case is for j > 0. The spectrum is also periodic in this case, i.e. it is the same at θ and θ + π, but the flow shifts the energy levels by one unit each time that we increase θ by π. In general, for fixed angular momentum and charge we have for any integer k ∈ Z. But even in that case the spectral flow has a monotonic character, i.e.
The inequalities between the Hydrogenoid and meta-Hydrogenoid energy levels (56) become the standard inequality of the Hydrogenoid levels E H n < E H n+1 in this case. The behaviour of the spectral flow recalls again the pumping mechanism of edge states in topological insulators [48,49].
The degeneracy between the bound energy levels with total angular momentum j and −j (for n > 0) at θ = 0 and θ = π is again broken for intermediate values of θ ∈ (0, π) as shown in Figure 3. The level with negative angular momentum −|j| has always lower energy than that of the corresponding level with positive angular momentum |j|.
Let us now analyze the supercritical charge regime with α 2 > j 2 . As anticipated, in this regime, the levels E H n and E h n do not belong to the spectrum of the Hamiltonian. In this case for any value of θ the energy spectrum E IV n (θ) contains an infinity number of bound states that accumulate near the mass gap continuum energy level E = m. In Figure (4) we plot the flow of some eigenvalues of the spectrum E IV n (θ) when parameter θ flows from 0 to π. Notice that along that flow one eigenvalue pops up from the Dirac sea continuum E < −m at a particular value of the parameter θ.
This is the only footprint of the instabilities pointed out in the supercritical regime, where the analytic expressions of Hydrogenoid and meta-Hydrogenoid energy levels become formally complex. Notice that the same phenomenon occurs in the subcritical regime α 2 < j 2 . The appearance of these instabilities is what inspired the Gribov mechanism of quark confinement in QCD [10]- [11] (see also [12]- [15]).

Spectral flow and impurity charges
In order to analyze the transition from the subcritical regime to the supercritical regime we fix a suitable value of the parameter θ for E II n (θ), E III n (θ) and E IV n (θ). By increasing the value of α we can follow the flow of each energy level from the subcritical regime to the critical regime in an adiabatic continuous way. Notice, however, that for each 0 < θ < π 4 there is a bound state in the subcritical regime that merges to the continuum for a special value of α < j , and conversely, there is an infinity of bound states emerging from the continuum spectrum for α |j| in the supercritical regime for any θ. In any case we have the following relations whenever θ = π 4 and θ = π 2 . This means that, for any fixed values of θ (θ = π 4 and θ = π 2 ), E II n (θ) and E IV n (θ ), converge to E III n (0) as α → α = |j|, pointing out the continuity of the flow of energy levels in the transition from the subcritical regime to the critical one (See Figure 5). In the exceptional cases we also have continuity in the path crossing the transition border  They are also highly dependent on the boundary conditions of the different self adjoint extensions. For simplicity we consider only the angular momentum j = 3 2 . For simplicity, only the flow of lowest bound states of the infinite tower is displayed for different values of α. The flow of the higher energy levels is in fact very similar. In the region 0 < α < √ 2, the operator is essentially self adjoint and the spectrum is that of an Hydrogenoid atom (black lines in Figure (5)) that begin at E = m for α = 0. On the border of the subcritical region (α = √ 2) the smallest level (red point in Figure (5)) is the ground state of the meta-Hydrogenoid spectrum, while the other two levels are doubly degenerated, because they include the n-level of the Hydrogenoid spectrum and the (n + 1)-level of the meta-Hydrogenoid spectrum. All the energy levels of the different self adjoint extensions of H start from one of these points.
At the critical coupling α = 3 2 , we also have some special energy levels, which to some extent, are attractors or repulsors of the other energy levels: the black points correspond to the double degenerate Hydrogenoid and meta-Hydrogenoid spectra of θ = 0, whereas the green points correspond to the bound states of the spectrum of H for θ = π/2. The alternating black and green points are, respectively, stable and unstable fixed points for the flow of energy levels. Each green point attracts only one energy level, that corresponds to the self adjoint extension with θ = π/4 which are on the green curve of Figure (5). These flow lines are isolated and act as repulsive barriers creating bifurcations of the flow. For 0 ≤ θ < π/4 the ground state merges into the continuum flowing to −∞ and all other n + 1 levels flow into the n black point, while for π/4 < θ < π all n levels flow towards the n black point (n ≥ 0). In the supercritical region, each green point is the starting point of only one level (green curve) which is associated to a particular selfadjoint extension θ = π/2. whereas black points are the starting points of bound state energy levels for all other boundary conditions θ = π/2. The green levels again are isolated and create a barrier for all the others. Notice that for any θ = π/2 there are energy levels which emerge from the continuum for large enough values of the charge α. In fact there is an infinity of them if we consider higher excited bound states. In Figure (5) we just displayed one of those levels emerging from the continuum for each boundary condition). The colors correspond to different choices of boundary conditions. In the subcritical regime θ = 0 (Hydrogen), θ = 0.005 π, θ = 0.03 π, θ = 0.1 π, θ = 0.2 π, θ = 0.25 π (isolated), θ = 0.3 π, θ = 0.5 π (meta-Hydrogen), θ = 0.9 π, θ = 0.99 π, θ = 0.999 π; and in the supercritical regime θ = 0, θ = 0.1 π, θ = 0.25 π, θ = 0.5 π (isolated), θ = 0.6 π, θ = 0.75 π, θ = 0.9 π.
In Figures (6) and (7) we show the instability of the isolated flow lines. The central flux lines correspond, respectively to θ = π 4 and θ = π 2 , while the others correspond to small perturbations of these lines, respectively θ = π 4 ± 0.001 and θ = π 2 ± 0.005. We can see how, when approaching to α 2 = j 2 , the perturbed curves follow the isolated lines flow but eventually they are attracted by two different eigenvalues of E III n (0).

Conclusions
In summary, the above analysis in terms of boundary conditions shows that in graphene we have infinite set of self adjoint Dirac operator for any α > 0 (for j = 1 2 ) which are parameterized by an ultraviolet scale Λ and an angle θ ∈ [0, π). In this sense the behaviour of impurities in graphene is different from that of Hydrogenoid atoms in QED. The lowest angular momentum states in graphene are always in the subcritical regime unlike in the 3D Hydrogen atom, which requires the introduction of appropriate boundary conditions that depend on a UV scale λ and a dimensionless parameter θ.
Unitarity is guaranteed for any value of the charge impurity α. Even more, the parameters introduced by the boundary conditions Λ, θ that renormalize the singular UV of the impurity induce remarkable observable effects. The dependence on the choice of boundary conditions at the singularity defines a flow of energy levels. The analysis of the flow of boundary levels displays interesting physical properties. Changes of the θ parameter which characterizes the self adjoint extension of the Hamiltonian can pump each Hydrogenoid level into the next one after a recursive loop in the parameter space recalling the pumping mechanism of topological insulators. A Berry phase can be also associated to this process. All energy levels in the Hydrogenoid spectrum, except the fundamental one, are degenerate, but the introduction of the parameter θ breaks down this degeneracy. Moreover, it is possible to change, by adiabatical variations of α, the energy levels from the subcritical to the supercritical regime in a continue way. Some bound states emerge (merge) from the continuum in this process. This is a consequence of the interesting properties of the RG flow for the subcritical and supercritical regime. Near the critical charge the energy levels are attracted by the points of the spectra of the Hamiltonian at the critical charge α 2 = j 2 and the particular value of the boundary conditions θ = 0. The attracting Hamiltonian corresponds to the Hydrogenoid atom spectrum at the critical charge. Only few levels remain isolated in a unstable way. This points out that the critical charge α 2 = j 2 of the Hydrogenoid case is not a singular case from the quantum physics viewpoint. The theory is well defined below and above this critical charge in the subcritical and supercritical regimes. The transition from the subcritical to the supercritical regime does not imply a critical change in the physical description of the system. However, the apparent stability of the vacuum pointed out by the careful analysis of the boundary conditions of the Hamiltonian can not hide that the physical behaviour of graphene is quite special in the supercritical phase. The fact that Hydrogenoid energy levels become complex in the supercritical regime implies the presence of resonances in the spectral density of the scattering matrix in the positron (hole) channel. These resonances are also the root of bound states levels which emerge from the continuum negative spectrum E < −m (see Figure  4).
In the supercritical regime there is an infinite number of quasi-bound states embedded in the lower continuum E < −m which are visible in the spectral density. If they are not filled when cross the supercritical value, some normal electrons will jump into these empty levels generating particle/hole pairs. The positive charges will move to infinite and disappear whereas the negative charges remain localized near the impurity giving rise to a screening of the impurity charge. We have assumed a positively charged impurity but due to the CP invariance of the theory a similar phenomenon occurs for negative charged impurities.
The phenomena described above are reminiscent of what happens in Quantum Eletrodynamics [5,6]. The main difference is that the value of the critical charge in graphene is α = j whereas in QED is Z = 137, which is very hard to realize in Nature. The screening phenomenon due to supercritical pair creation has been recently observed in graphene [19,18] and in QED a similar phenomenon might be also observed in the heavy ions collisions (see [50] for an updated review). There is another remarkable difference between the two theories. In graphene for any value of α > 0 the system is in a subcritical regime at least in the lowest angular momentum sector (j = 1 2 ) which requires always the choice of an extra parameter to fix the boundary condition at the origin. However, in QED for Z < 137, e.g. for the Hydrogen atom the Hamiltonian is essentially selfadjoint in lowest angular momentum sector. Thus there is no need to fix the boundary condition at the origin. In particular a potential δ like perturbation has no effect in the spectrum. In particular this means that the relativistic interpretation of the Lamb effect cannot be understood in pure relativistic quantum mechanics and requires a full field theoretical analysis, unlike in the non-relativistic quantum mechanics approach.
The analysis of the energy spectrum in the gapless semi-metal regime of graphene can be carried out in a similar way. The boundary conditions are exactly the same as in the massive case and, thus, the different physical regimes are the same. However, there is a fundamental difference, there are no electronic bound states, because the exponential decay at infinity disappears when m → 0, although there are some special points of the continuous spectrum that correspond to resonances which can be observed in scattering processes [18,19].