{'fraction_non_alphanumeric': 0.05195651529962381, 'fraction_numerical': 0.017387955804378092, 'mean_word_length': 3.924872355944566, 'pattern_counts': {'":': 0, '<': 18, '<?xml version=': 0, '>': 26, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Bound and resonance electron states in the monolayer graphene with the short-range impurities

Bound and resonance electronic states in impure graphene are studied. Short-range perturbations for defects and impurities of the types"local chemical potential"and"local gap"are taken into account. Zero gap and non-zero gap kinds of graphene are considered. A qualitative analysis of the electronic spectrum general features and numeric calculations are presented. A dependence of the resonance widths on the gap value, angular momentum values and the perturbation amplitudes is investigated.


I. INTRODUCTION
The Dirac equation is one of keystones of the relativistic field theory. However, it is an important model in the non-relativistic solid state theory as well. Superconductors with d−pairing [1], the Cohen-Blount two-band model of narrow-gap semiconductors [2], [3], electronic spectrum of the carbon tubes form an incomplete list of the nonrelativistic applications of this equation. During the last two years extremely much attention was payed to the problem of the electronic spectrum of graphene (see a review [4]). Two-dimensional structure of it and a presence of the cone points in the electronic spectrum make actual a comprehensive study of the external fields effect on the spectrum and other characteristics of the electronic states described by the Dirac equation in the 2+1 space-time. We consider in this work the resonance and bound states of the 2+1 Dirac equation due to the short-range perturbation using some approaches developed in Ref. [5] .We do not take into account the inter-valley transitions. Particular attention to this case stems from the effectiveness of short-range scatterers in contrast to the long-range ones: an effect of the latter is suppressed by the Klein paradox [6]. Our work takes into account the obvious fact that the Kohn-Luttinger matrix elements of the short-range perturbation calculated on the upper and lower band wave functions are not equal in a general case. This means that not only the potential but the mass perturbation can be present in the perturbed Dirac equation. Resonance states in the zero-mass Dirac model of graphene were considered for distinct perturbations in [7], [8], [9] .
The paper is organized as follows. In the first section we formulate the problem in terms of the 2+1 Dirac equation with the delta function perturbation taking into account a local change of the potential and mass (gap) due to the crystal defect. An exact solution of this problem leads us to the characteristic equation determining the bound and resonance states in graphene. In the second section we analyze the characteristic equation for the zero-gap graphene case. Both local and global behavior of resonance states for various magnitudes of the perturbations are investigated. Analytic and numeric results are presented. Just analytic approach allowed us to investigate the case of exponentially narrow resonance that could be difficult to do numerically. In the third section we consider the bound and resonance states in graphene with non-zero mass (gap). The exponential approach of the bound state energy levels to the band edge was studied analytically; numerical analysis is presented as well.

II. CHARACTERISTIC EQUATION
The Dirac equation describing electronic states in graphene reads [4] where v F is the Fermi velocity of the band electrons, γ µ are the Dirac matrices γ 0 = σ 3 , γ 1 = σ 1 , γ 2 = iσ 2 , σ i are the Pauli matrices, 2mv F 2 = E g is the electronic bandgap, ψ (r) is the two-component spinor. The electronic gap can appear in the graphene monatomic film lying on the substrate because of the sublattices mutual shift [10]; the spin-orbit interaction can be the reason too. The spinor structure takes into account the two-sublattice structure of graphene. δm (r) and V (r) are the local perturbations of the mass (gap) and the chemical potential. A local mass perturbation can be induced by defects in a graphene film or in the substrate [10]. We consider here the delta function model of the perturbation: δm (r) = −bδ(r − r 0 ), V (r) = −aδ(r − r 0 ), (2) where r and r 0 are respectively the polar coordinate radius and the perturbation radius. Such short-range perturbation was used in the (3+1)-Dirac problem for narrow-gap and zero-gap semiconductors in [3]. The perturbation matrix elements are related to the a, b parameters as follows The delta function perturbation is the simplest solvable short-range model. Finite radius r 0 plays a role of the regulator and is necessary in order to exclude deep states of the atomic energy scale. The finite perturbation radius r 0 leads to the quasi-momentum space form-factor proportional to the Bessel function that justifies our neglect of transitions between the points K and K ′ Ref. [3]. The two-dimensional Dirac problem with the scalar shortrange perturbation Eq. (2) (but without the mass perturbation) was considered in Ref. [11]. The obtained there characteristic equation for the discrete energy spectrum contains one mistake. We corrected it in our previous work Ref. [12] and took into account the mass perturbation δm (r) . Here we present more detailed analysis of the electronic bound states and new results on resonance states in the zero-gap and non-zero-gap graphene. Particular attention is paid to the effect of relative intensities of two perturbations V 1 and V 2 (a and b).
Let us present the two-component spinor in the form where j is the pseudospin quantum number; j = ±1/2, ±3/2, . . .. In opposite to the relativistic theory, this quantum number has nothing to do with the real spin and indicates a degeneracy in the biconic Dirac point. The upper f j (r) and lower g j (r) components of the spinor satisfy the equations set These equations have a symmetry: Let us introduce the function ϕ j (r) ≡ f j /g j . It satisfies the equation: Integrating in the vicinity of r = r 0 lim ǫ→0 ϕ j (r0+δ) we obtain the matching condition where ϕ − j ≡ ϕ j (r 0 − δ) , ϕ + j ≡ ϕ j (r 0 + δ) , δ −→ 0 . The upper and lower component matching conditions resulting from Eq. (11) read where is the orthogonal rotation matrix. It transmutes into the orthogonal boost matrix for b 2 > a 2 The general solution can be found solving the second-order equation obtained by excluding one of the spinor components from the equation set Eq. (6), Eq. (7) in the domains 0 < r < r 0 and r > r 0 : This equation is related to the Bessel one. We assume at first E to be real and satisfying the inequality E 2 < m 2 .
Then the general solution of Eq. (15) reads where κ = √ m 2 − E 2 is the principal value of the root; I ν (z) and K ν (z) are the modified Bessel functions. The constant C 2 = 0 in the domain 0 < r < r 0 , while C 1 = 0 in the domain r > r 0 . It is useful to introduce the following notations: Expressing the g j -component from Eq. (7), we can write Substituting the expressions Eq.(18), Eq.(19) into the matching condition Eq. (11), we obtain the characteristic equation for the bound state energy levels (for E 2 − m 2 < 0): where T (a, b) .is determined as follows: We presented above an analysis for the case of a 2 > b 2 ; the case of b 2 > a 2 can be considered similarly.
This equation turns to the characteristic one obtained in [11], for b = 0 apart from the mistakenly omitted terms in the right hand side of Eq. (20). This characteristic equation is unambiguously determined for the bound states with energy levels lying in the real axis segment [−m, m] .
We write this equation in another form making the symmetry Eq. (8) manifest: This equation was derived for a study of bound states situated in the gap that will be done in section IV. An analytic continuation of this equation from the energy real axis segment (−m, m) onto the energy complex plane is necessary in order to study resonance states. Using the obvious relations we obtain the characteristic equation in the form where We have introduced the notation: The cases of a 2 = b 2 (coordinate angles bisectrices in the (a, b)-plane) are degenerate: the energy imaginary part vanishes in this limit since according to Eq. (4), one of the matrix elements V 1(2) equals zero. For instance, the effective second-order equation for the upper spinor component takes the simple form: This equation obviously does not contain a resonance state because of absence of terms of the type characteristic for the relativistic resonances V 2 , dV /dr etc [14].

III. ZERO-GAP GRAPHENE
Let us consider Eq. (25) in the zero-gap case m = 0. The symmetry Eq. (8) allows us to restrict the analysis by the right energy half-plane. If we consider explicitly only the right energy half-plane (α = 1), we should take into account both positive and negative angular momentum quantum numbers j = ±1/2, ±3/2, . . .and both positive and negative potential amplitude a values in order to obtain a complete picture, while considering both α = 1 and α = 2 cases it was enough to take into account only positive j. The characteristic equation can be essentially simplified in the zero-gap limit m = 0 : Notice that this equation has a root E = 0. Vanishing of the energy imaginary part stems obviously from vanishing of the free Dirac density of states g (E) = T r Im G r (E) at E = 0. Here G r (E) is the free retard Green function. Complex roots of the characteristic equation will be interpreted here as resonance states. Now we will study a distribution of the complex roots of Eq. (28), i. e. resonances.
It is convenient to divide the plane (a, b) into four regions separated by the coordinate angles bisectrices a 2 = b 2 (see Fig. 1). We numerate these regions from I to IV.
A. Resonances: a 2 > b 2 Let us consider at first the region I of the (a, b)-plane (see Fig. 1). It is useful to introduce there the hyperbolic variables: Then Eq. (28) takes the form: Periodicity of tanρ 1 allows us to restrict the analysis here by the segment 0 ≤ ρ 1 ≤ π. Consideration of only right half-plane of energies Re E > 0 and j = ±1/2, ±3/2, . . .will be suffice due to the symmetry Eq. (8).
For a fixed ρ 1 , we can consider the asymptotic behavior at |ψ 1 | → ∞. Eq. (30) asymptotically transforms into the equations Since the Hankel functions H j∓1/2 (Er 0 ) have no roots in the right halfplane [13], solutions E (j∓1/2) n of Eq. (31), Eq. (32) read with λ (j∓1/2) n being the roots of the Bessel functions J j∓1/2 (z) , which are real and tend to infinity. Thus we conclude that solutions of Eq. (30) are distributed on the arcuated curves beared on the real axis in the points E (j∓1/2) n (see Fig. 2) Notice that the case of |ψ 1 | → ∞ corresponds to a 2 ⇆ b 2 (see section II, Eq. (27). When we move along the non-singular hyperbolae ρ 1 ∈ (0, π/2) from ψ 1 = −∞ to ψ 1 = +∞, the n-th arc-like curves are circumscribed in the complex plane of E from the lying on the real axis point E When ρ 1 approaches π/2 (singular hyperbola), the "arc" height tends to infinity. Moving along the arc-like curve takes place anti-clock-wise, when ψ 1 varies from −∞ to +∞. A transition between the points E (j±1/2) n corresponds to swapping of the upper and lower components of the wave function spinor amplitudes. However, crossing the singular hyperbolae π/2 leads to an essential rearrangement of the arc structure. Notably, the function tan ρ 1 changes its sign and alternative combinations of nearest neighbor Bessel's function roots are connected by arcs in order to satisfy the characteristic equation. Thus, the resonance arc-like curve begins in the point E (j+1/2) n and moves clock-wise into the point E , when π/2 < ρ 1 < π. A schematic picture of such transitions of arcs is presented in Fig. 3. Results of numerical calculations for ρ 1 ∈ (0, π) presented in Fig. 4 for ρ 1 ∈ (0, π/2) and in Fig. 5 for ρ 1 ∈ (π/2, π) . They are in good qualitative agreement with our analytic consideration. A periodicity with the period ∆ρ 1 = π takes place.
Motion along rays in the region I of the (a, b)-plane maps into closed resonance trajectories in the energy complex plane in the case of the ray lying in the upper half-plane. In the case of the ray lying in the lower half-plane we obtain repeating jumps from one resonance to the next one instead of the closed curves (see Fig. 6)..When the ray approaches the coordinate angle bisectrix, the closed curves diameter decreases so that the curves shrink into a point.
Consider now the left quarter (region III in Fig. 1): Then the characteristic equation reads The asymptotic form of the equation at |ψ 3 | −→ ∞ is following When ψ 3 varies from −∞ to +∞ circumscribing the hyperbola in the region III, we have motion from the roots of J j+1/2 (r 0 E) to the roots of J j−1/2 (r 0 E) in the energy complex plane (the same arc-form curves circumscribed anti-clock-wise or clock-wise respectively for ρ 3 ∈ (0, π/2) or (π/2, π)). Numeric calculations are similar to the case of the region I; they confirm conclusions of our analysis. An increase of the potential radius r 0 shifts the roots in the direction of the coordinate origin.
Now we consider the region II in the (a, b)-plane (see Fig. 1). We introduce the hyperbolic variables here differently: The characteristic equation takes the form One can see that in the limit of a ⇆ b we have the equation Thus, the roots circumscribe arc-form curves clock-wise from the point λ (j−1/2) n /r 0 to λ (j+1/2) n /r 0 (see Fig. 7), while ψ 2 varies from +∞ to −∞ (motion along the upper hyperbolae in Fig. 1 anti-clock-wise).
Let us learn a limiting form of the arcs, when ρ 2 −→ ∞. As tanh ρ 2 −→ 1, Eq. (38) takes the form where E jn = E jn − iΓ jn /2 is the complex root of the characteristic equation; we used the standard notions for the real and imaginary parts of the resonance energy.
Notice that roots of the second factor are not physical and must be ignored. Eq. (39) can be written in the form Since J j (z) is real on the real axis, we have J j (z) * = J j (z * ) and, therefore, jn + iΓ jn /2 is a root of Eq. (40) as well. Since the ends of the arc are simple zeros of the Bessel function, these two arcs coincide and are flat in this limit. Eq. (40) takes the form We see that the energy imaginary part is getting extremely small and the arcs approach the real axis segments to observe a presence of the imaginary part numerically. However, this process can be followed analytically. In order to investigate the limiting behavior of the energy imaginary part, we expand Eq. (38) for ρ 2 −→ ∞ taking account of the asymptotics leading terms. In particular, the expansion tanh x ∼ = 1 − exp (−2x) + . . . for x ≫ 1 was used. Taking the real part of Eq. (38), we obtain after some algebra Expanding this formula for small Γ jn we obtain in the leading order Here E jn is determined from Eq. (41). It is seen from Eq. (42) that the resonance width Γ jn decreases exponentially at large ρ 2 .
Returning to the initial notions we find: The obtained formula indicates that these resonances are extremely sharp for the case of small |j| and, therefore, they can give a strong resonance scattering of electrons. Using the Bessel functions J ν , H (1) ν asymptotics for large ν [13], we obtain the following estimate for the resonance width asymptotic at large |j| : This function catastrophically increases beginning from |j| ∼ 4. It means that the imaginary part of the energy is large for any realistic ρ 2 so that solutions with j > 4 are not really sharp resonances and they will not contribute into the electron resonance scattering.
Results of numeric calculation for resonances in the case of b 2 > a 2 are presented in Figs. 7 and 8. Approaching of arc-like resonance curves to the real axis segments, when ρ 2 increases, is shown in Fig. 7. Fig. 8 shows the resonance trajectory corresonding to a motion along the ray in the (a, b)-plane region II. When the radius increases, the trajectory approaches a point on the real axis.
The region IV (lower quarter) can be analyzed similarly using the hyperbolic variables: An analisis of this region could be done similarly to the region II, but we omit it here since no qualitative distinction is observed in this case. Numerical results for the region IV are omitted by us as well because of this similarity.
So we have discussed a local behavior of the resonance trajectories within the chosen quarter in the (a, b)-plane. Now we will consider the global behavior, when all quarters are passed along hyperbolae with the fixed hyperbolic radius ρ. This process creates pairs of arc-form curves bearing on two common basing points lying on the real axis. These figures of eight are arranged at the E complex plane symmetrically relative to the imaginary axis in the case of m = 0. Resonances move anti-clock-wise from the root λ (j+1/2) n /r 0 to the root λ (j−1/2) n /r 0 circumscribing the arc, when φ ∈ [−π/4, π/4] When φ ∈ [π/4, 3π/4] , a different arc appears. It connects the same points λ (j+1/2) n and λ (j−1/2) n , but moving in the opposite direction (clock-wise). We obtain in result the figures-of-eight in the E-complex plane. This picture is repeated with the period ∆φ = π.
Such excursion around all quarters of the (a, b)-plane can be carried out circumscribing the circles in the (a, b)-plane One quarter after another is passed in this case too, but now the figures-of-eight do not touch the real energy axis anymore (see Fig. 9), i. e. all resonances belonging to this set have rather large widths. Numeric analysis confirms our prediction on the mapping of the circular excursion in the (a, b) plane onto the energy plane. When the perturbation intensity ρ is increased, the picture is getting more complicated: all the figures-of-eight are associating into a single whole (see Fig. 9 for large ρ). Notice that the resonance trajectory self-crossing points are fixed points relative to ρ changes. !They correspond to resonance set, which do not change at motion along the ray a = b.

IV. GAPPED GRAPHENE
While the pristine graphene is gapless, a contact with the substrate can induce some narrow gap [10]. That is why we study here bound and resonance states in the gapped graphene as well.
Let us consider the characteristic equation for the case of m = 0. We introduce the dimensionless parameters: where r c = 1/m is the "Compton radius" (ℏ = 1, v F = 1). The characteristic equation for bound states Eq. (22) takes the form: The perturbation radius r 0 value is of order of the lattice spacing d, while the "Compton" radius r c is rather large: r c >> d : ℏ/ (m 0 v F ) ∼ 10 −8 cm even for the bare electron mass m 0 ; the real effective mass at the band extrema can be estimated as 10 −1 − 10 −2 of m 0 . Therefore, the actual value of γ is less than unity.
We consider here some general properties of the bound states electronic spectrum resulting from the characteristic equation (48) and the effect of non-vanishing gap on the resonance states. Results of the numerical solution will be presented too.
A. Bound states: a 2 > b 2 Let us assume at first j = 1/2. We consider Eq. (48) for the energy situated in the gap ǫ 2 < 1. We use the Bessel functions limiting forms for the states lying near the gap edges [13] I ν (z) ∼ (z/2) in order to transform the characteristic equation (48). Here Γ (ν) is the gamma function. We obtain a simple relation for small z.
Formulae like Eq. (50) determine a map from the plane (a, b) onto the complex energy plane; in the case of bound states the map is carried out onto the real energy axis segment −1 < ǫ < 1.Now we consider the asymptotic energy level behavior near the upper edge of the gap ǫ −→ 1. Eq. (50) takes the asymptotic form: Therefore, we can write It is clear that ǫ −→ 1 if a+b −→ +. 0. This result conforms the well known general property of the two-dimensional quantum systems: a threshold for creation of the bound state is absent for j = 1/2. Formula (52) can be rewritten in the hyperbolic variables in the region I (see Eq. (29) as follows It is seen from the last formula that ǫ −→ 1, when ρ 1 −→ πn + 0 (where n = 1, 2 . . .), i. e. the ǫ = 1 "fronts" approach the hyperbolae ρ 1 = √ a 2 − b 2 −→ πn + 0. Let us study now a behavior of the bound state trajectory near the lower band edge. Using Eq. (50) we conclude that when ǫ −→ −1, the following relation holds: Other domains with n = 1, 2, . . . , lie between approaching the hyperbolae ρ 1 = πn − 0 curves (ǫ = 1 "front", |ψ 1 | −→ ∞) and the ǫ = −1 "front", approaching the hyperbolae ρ 1 = π/2 + πn − 0 at ψ 1 −→ +∞ and the hyperbolae ρ 1 = πn + 0 at ψ 1 −→ −∞. These lines boarder the countable set of domains in the (a, b)-plane, where the bound states exist. When we move along the lying in the region I ray beginning in the coordinate origin in the direction of the ρ 2 1 = a 2 − b 2 increase, we cross this set of domains so that every time the bound states trajectory starts at ǫ = 1 and terminates at ǫ = −1.
A dependence of the bound state energy on the radius ρ 1 , at the motion along the rays in the (a, b)-plane (see inset to Fig. 10) was investigated numerically in the case of a 2 > b 2 . The motion between the coordinate origin and the first ǫ = −1 "front" maps into the mildly sloping transition between the ǫ = 1 and ǫ = −1 edges of the gap (solid line) and the sharp sloping transition respectively for the rays in the upper and lower half-planes of the (a, b)-plane (first pair of the curves). The second and so on pairs of curves show similar transitions between the domain boundaries. Now we consider the region III. Let us return to Eq. (50) that is valid near the band edges and everywhere inside the gap if γ << 1. Using the hyperbolic variables Eq. (35) we can write for the bound states near the upper edge ǫ = +1 : We obtain for the energy near the upper band gap edge: Thus, the ǫ = 1 "front" approaches the bisectrix a + b −→ 0 exponentially. We introduce here the variables Then Eq. (58) takes the form We obtain for the energies near the upper band edge ǫ −→ 1: Thus, the ǫ = 1 "front" approaches the bisectrix a + b = 0 at ρ 2 −→ 0. There is no "front" at the bisectrix a − b = 0 as it was shown above. Eq. (60) can be rewritten for ρ 2 −→ ∞ : Therefore, Eq. (63) is satisfied if ǫ −→ −1 and ψ 2 −→ +∞. Therefore, the ǫ = −1 "front" approaches the bisectrix a − b = 0, if ρ 2 −→ ∞.This "front" lies in the region I (see above). Thus, a motion in the region II of the plane (a, b) from the bisectrix a + b = 0 to the bisectrix a − b = 0 corresponds to a motion of the energy eigenvalue from the upper band edge ǫ = 1 to the lower band edge ǫ = −1, when ρ 2 >> 1.
A dependence of the bound state energy on the radius ρ 2 , at the motion along the rays in the (a, b)-plane (see inset to Fig. 11) was investigated numerically in the case of a 2 < b 2 . The motion along the ray maps into the monotone transition of the energy level from the band gap edge ǫ = 1 to the constant situated in the gap. Motion along the ray in the right half-plane maps into the faster energy dependence and with lower asymptotic energy value, than in the case of the ray in the left half-plane. When the ray approaches the bisectrix a − b = 0, the asymptotic value tends to ǫ = −1.
Let us consider now solutions of the characteristic equation near the band edges for the parameters a, b in the region IV. We introduce the hyperbolic variables (see Eq. (45)).
Then Eq. (58) can be approximately written in the form: Let us transform Eq. (64) in the limit of ρ 4 −→ ∞ : The "front" ǫ −→ 1 equation takes the form Thus, the "front" ǫ −→ 1 approaches the bisectrix b + a = 0 asymptotically, when ψ 4 −→ −∞. Therefore, a motion from the line Eq. (66) to the bisectrix b + a = 0 in the (a, b) plane is mapped to a motion from the lower edge of the gap to the upper one when ρ 4 >> 1.
If γ increases, we conclude from Eq. (66) that the second front 1 + ǫ −→ 0 moves to the bisectrix b − a = 0, i. e. an increase of the mass widens the domain of the region IV, where the eigenvalues exist; it fills all region IV in the extreme limit.of large γ. Notice that this case is not realistic for an impurity problem, but can be important for the quantum dot case [9].
Thus, our analysis shows that there is a countable set of eigenvalue domains in the case a 2 > b 2 . The boundaries of these domains approach the hyperbolae described above. In the case of b 2 > a 2 , we have a saturation instead of periodicity, when radius increases.

C. Bound states for higher angular momenta
Let us consider now the eigenvalue spectrum for the angular momentum j = 3/2,for instance, when the parameters a, b lie in the regions I and II.. Using the Bessel functions expansion for small arguments (49), Eq. (48) can be written as follows: In the region I, this equation takes the form The "fronts" ǫ −→ 1 have the form We have asymptotically: The "fronts" ǫ −→ −1 have the form tan ρ 1 = (6/γ) exp (ψ 1 ) .
Therefore, a motion from the point a 1 to the point a 2 maps into a motion from ǫ = 1 to ǫ = −1. More than that, having in mind a continuous dependence of solutions on the parameters a and b the "fronts" Eq. (71), Eq. (73) obtained from the asymptotic equation (70) cross one another. This prediction is confirmed by the numerical calculations (see Fig. 13). Exclusively fast dependence of the energy eigenvalues on the parameters a and b will be observed in the vicinity of the crossing point. Let us consider now the region II. We have here from Eq. (69): The "front" ǫ −→ 1 is given by the equation Since 0 ≤ tanh ρ 2 ≤ 1 (ρ 2 ≥ 0) , we have The domain, where the energy eigenvalues exist in the case of γ = 2, is bounded by the positive half-axis b and the first quadrant bisectrix. When γ < 2 and decreases, the "front" ǫ = 1 turns about the coordinate origin from the ordinate axis approaching the first quadrant bisectrix. The eigenvalues domain is narrowing in result. When γ > 2 and increases, this boundary turns anti-clock-wise approaching the second quadrant bisectrix. When ψ 2 −→ ∞, Eq. (77) can be satisfied only if We have for ρ 2 −→ ∞ (higher hyperbolae): When ψ 2 is fixed, the asymptotic magnitude of the eigenvalue ǫ as at ρ 2 −→ ∞ can be found from the quadratic equation The domain of eigenvalues in the region IV is narrower in the case of j = 3/2, than for j = 1/2 similarly to the region II. Generically, an increase of j at fixed γ leads to narrowing of the eigenvalues existence domains.
The case of higher angular momenta j = 5/2, . . . can be analyzed similarly. The energy discrete eigenvalues exist for all j, but in opposite to the case of j = 1/2, there is a threshold for the bound state appearing at the upper edge of the gap. Existence of bound states for all higher angular momenta is a specific feature of the delta function perturbation.
The states with j < 0 can be easily obtained with a use of the symmetry transformation Eq. (8).

D. Diagram of the bound states
General results of mathematical and numerical analysis for a distribution of bound states for j = 1/2 and j = 3/2 are presented respectively in Figs. 12, 13. Shaded domains represent basines of the bound states existence. The thick lines correspond to the upper edge of the band gap ǫ = +1, while the thin ones represent the lower edge ǫ = −1.
In the case of j = 1/2, domain boundaries approach asymptotically the hyperbolae. One crossing point of boundaries can be seen. Therefore, we have one twisted domain in this case. Notice that this crossing looks so in a projection onto the (a, b)-plane. A three-dimensional picture of this feature is shown.in Fig. 13: it is seen that the "fronts" ǫ = ±1 lie at different "floors" so that the crossing curves in the plane curves appear to be skew ones.
Notice that our diagrams of states are obtained assuming single-valued dependence of the bound state energy on the parameters a and b. This is true everywhere in the (a, b)-plane except very small islands near the domain boundaries (see Fig. 10 of this paper and our work [12], where the twisting point neighborhood was analyzed numerically).
In the case of j = 3/2,the bound state existence domains are getting narrower, a number of the boundaries crossing points increases. The diagram of states for higher j is given us in Ref. [12].

E. Resonances: m > 0
It is necessary to return now to Eq. (25) written in terms of non-modified Bessel's functions. We consider here an effect of non-zero bandgap upon the resonance states. It is convenient to re-write Eq. (25) using the Bessel function explicitly (see eq (26)) and choosing α = 1. We analyze here resonances for the parameters values lying in the region I as an example.
Introducing the hyperbolic coordinates according to Eq. (29) and using the dimensionless parametrization Eq. (47) , we can write Re E > 0, j = ±1/2, ±3/2, . . . , where ς = γ √ ǫ 2 − 1. If ψ 1 −→ ±∞ (asymptotical approach to the bisectrix), we have for the roots since the roots of the equation H j±1/2 (ς) = 0 are extraneous ones. Solving Eq. (84), we obtain the Bessel function roots λ (j∓1/2) n . Therefore, This formula can be written in the form: We see that appearing of non-zero mass shifts "arc's" bearing points relative to the positions, determined by Eqs.. (33), (34), in the direction of the energy increase (see Fig. 14). This shift is essential only for small n since λ (j∓1/2) n increases fast with n increasing. When m = 0, a distribution of signs of Eq. (83) right-hand-side for various m values at fixed ρ 1 is complicated. Arcs, in result, can be formed both clock-wise and anti-clock depending on the m value.
We have analyzed the resonances for the region I; the dependence of the resonances distribution on the parameter γ = r 0 m in other regions can be considered similarly.
We studied in the previous section a motion of resonances along the figure-of-eight curves in the case of m = 0, when the parameters a, b are varied along the circles Eq.(46) with various ρ values. In the non-zero gap case, we have found that a large enough mass destructs the figure-of-eight structure for ρ fixed; these figures associate into a whole complicated aggregate (see Fig. 15).

V. CONCLUSION
In conclusion, we considered the bound and resonance electronic states for the two-dimensional Dirac equation with the short-range perturbation. The short-range perturbation is approximated by the delta function δ (r − r 0 ) with different amplitudes in the upper and lower electronic bands. We have in result local perturbations of the potential-like and mass-like types.respectively with amplitudes a and b.
We derived the characteristic equations for the bound and resonance states. The characteristic equation is presented in forms convenient both for an analysis of the bound and resonance states. Energy levels behavior in dependence on the perturbation amplitudes was investigated both analytically and numerically. A general picture of the electronic spectrum dependence on the potential, and mass local perturbation is considered for various mass and angular momentum values; the results are presented as diagrams in the (a, b)-plane. Absence of a threshold for forming of the bound state in the vicinity of the upper bandgap edge is obtained for the lower angular momenta |j| = 1/2 which is in accord with the general principles of quantum mechanics for the dimension 2+1.
A qualitative distinction of a 2 > b 2 and a 2 > b 2 cases is shown. A countable set of eigenvalues existence domains in the (a, b)-plane is present in the former case, while monotonic approaching to the asymptote takes place, when b 2 − a 2 −→ ∞. in the latter one. Positions of asymptotes are determined by the parameter tanh −1 (b/a) . Higher angular momenta are investigated as well. Twisted eigenvalue domains with crossing boundaries are found. Domains of eigenvalues are narrowing with the angular momentum increase and at the mass value fixed. A zero-energy solution with vanishing imaginary part exists in the zero-gap case.This solution has to be considered as a limit at m −→ 0.
Behavior of resonances depends essentially on the sign of a 2 − b 2 . When a 2 > b 2 , resonance trajectories in the energy complex plane have a form of closed curves circumscribed periodically at moving from the coordinate origin to infinity in the (a, b)-plane. Motion along the hyperbolae √ a 2 − b 2 = const maps onto the countable set of arc-like trajectories in the energy complex plane. These arcs bear on the real axis; positions of the bearing points are exactly determined by the Bessel function roots. When b 2 − a 2 > 0, an increase of b 2 − a 2 forces the arcs to get flatter so that if b 2 − a 2 −→ ∞, the arcs asymptotically nestle up to the real axis. The resonance width is exponentially small in this case. For higher angular momenta j, the rate of this decrease is getting small due to the factor |j| |j| . Therefore, higher spherical harmonics do not contribute into the resonance scattering.
Two kinds of the non-zero mass effect on the resonances behavior were found. Firstly, positions of arc-bearing points shift. This effect is particularly essential for lower roots. Secondly, eight-like figures resulted from mapping of the circular motion in the (a, b)-plane onto the energy plane transform into large aggregates.
The carried out analysis allowed us to describe both qualitatively and quantitatively a distribution of the bound and resonance states in the energy complex plane for any chosen distribution of the perturbation amplitudes a and b and for various magnitudes of the mass and momentum. The obtained results can be useful for understanding of the graphene electronic properties.    FIG. 3: Recommuting of the arc system in the region I, when ρ 1 crosses the singular hyperbolae ρ 1 = π/2 + πn, n = 0, 1 . . . (a) πn < ρ 1 < π/2 + πn.