Critical behavior for point monopole and dipole electric impurities in uniformily and uniaxially strained graphene

We revisit the problem of bound states in graphene under the influence of point electric monopole and dipole impurity potentials extended to the case in which the membrane of this material is uniformly and uniaxially strained, which leads to a redefinition of the charge and dipole moment, respectively. By considering an anisotropic Fermi velocity, we analytically solve the resulting Dirac equation for each potential. We observe that the effect of the anisotropy is to promote or inhibit the critical behavior known to occur for each kind of impurity, depending on the direction along which strain is applied: both the atomic collapse, for the monopole impurity, and the emergence of cascades of infinitely many bound states with a universal Efimov-like scaling, for the dipole impurity, are phenomena that occur under less or more restrictive conditions due to strain.


I. INTRODUCTION
Ever since the isolation of graphene samples in 2004 1 , this material has played a leading role in exploring analogs and similarities of phenomena occurring in different branches of physics. In particular, the low-energy quasi-particle excitations in this material, which behave as ultra-relativistic Dirac fermions in two space dimensions, have permitted to establish long standing predictions in quantum electrodynamics (QED) but in a condensed matter/solid state realm. Klein tunneling 2-4 , Zitterbewegung 5,6 and atomic collapse [7][8][9][10] are examples of phenomena that have been experimentally observed and theoretically predicted to occur in graphene due to the pseudo-relativistic character of its charge carriers near the Dirac points. In QED, the phenomenon of atomic collapse refers to the prediction that an electron would dive into the nucleus of a super-heavy element emitting a positron when the charge of the nucleus exceeds certain critical threshold, Z c 11,12 . The physical explanation supporting this unobserved prediction is that intense Coulomb field causes the electron wave function component to fall into the center of the nucleus (with a finite lifetime 13 ), and, at the same time, the positron component escapes to infinity. Such supercritical quantum states depict the semiclassical idea of an electron spiralling inward toward the nucleus as the corresponding positron spirals outward away from it, propagating to infinity. For a punctual nucleus, the threshold is calculated to be Z c = 137, whereas for extended nucleus, Z c 170 [11][12][13] . Up to date, this phenomenon has not been observed in nuclear or heavy-ion experiments nor in any other physical context.
In graphene, it has been observed by several groups that the low-energy Dirac Hamiltonian becomes non self-adjoint when the coupling of the long-range static Coulomb interaction for β = Zα/κ > β c 1/2 [7][8][9][10]14,15 . Here, α is the fine structure constant, κ the dielectric constant and Z the charge of an artificial nucleus. In this material, the critical charge Z c of artificial nuclei (Coulomb impurity) is strongly reduced [7][8][9]15,16 . This theoretical prediction has been experimentally confirmed 17 by using the tip of a scanning tunnelling microscope (STM) to create clusters of charged calcium dimmers, creating a supercritically charged Coulomb potential from subcritical charges (as in heavy-ion collisions). In these experiments, a systematic electron-hole asymmetry owing to the positive charge of the dimmers was observed, along with the appearance of a resonance (a quasi-bound state) increasing and shifting toward lower energies as the number of dimmers increased from one to five. For five dimmers, the resonance shifted below the Dirac point, representing the atomic collapse of the electron in graphene.
As a natural generalization, the problem of graphene under the influence of an electric dipole impurity has recently been addressed by several groups [18][19][20][21][22] . For a finite-size electric dipole impurity, it has been observed that when each of the charges of the dipole exceeds a critical value, the wavefunction of the highest energy occupied bound state changes its localization from the negatively charged impurity to the positive one as the distance between the impurities changes. This phenomenon has been dubbed as the migration of the electron wave function and comprises a generalization of the atomic collapse phenemonon by a single charged imputity to the scenario where both electrons and holes are spontaneously created from the vacuum in bound states with the two charges of the dipole partially screening them. Nevertheless, in the point electric dipole limit, which occurs in the limit when the distance between the charges of the dipole van-ishes with fixed electric dipole moment, a reminiscent of the critical behavior is observed: there exists at least one cascade of infinitely many bound states for an arbitrary value of the dipole strength. Additional cascades develop when the dipole strength achieves critical values according to a universal Efimov-like scaling [23][24][25] . There is no collapse of the electron wave function to the continuum because even though the electron and positron energy levels tend to come close to each other as the dipole moment increases, eventually they separate apart, thus preventing collapse.
In this article we revisit these critical problems in uniformly and uniaxially strained graphene. We consider the case in which the sample is first strained and then the impurities are later located in the plane. Mechanical deformations of the graphene membranes are known to modify the electronic properties of the material [26][27][28][29][30][31][32] . As a matter of fact, straintronics 33 has emerged as the field of research studying the modifications of graphene properties due to strain. For the simplest case we address in this article, it is well known that a uniform and uniaxial strain still induces an anisotropy in the Fermi velocity that, nevertheless, does not generate pseudomagnetic fields 31 . Thus, by solving the corresponding anisotropic Dirac equation under the influence of monopole and dipole electric impurity potentials, we explore the role of the anisotropy on the critical thresholds for these problems. For this purpose, we have organized the remaining of this article as follows: In the next Section we solve the anisotropic Dirac equation for the case of a monopole electric impurity immersed in a uniform magnetic field 34 and explore the behavior of the bound state energies. We regularize the Coulomb potential at the origin replacing it with a potential well extended up to a finite "elliptic radius" (this point will become clear in the discussion below) r 0 and of strength V 0 . We explore the effect of the anisotropy parameter in the collapse of bound states at different angular momenta for several regimes. In Section 3 we address the problem of a point electric dipole in uniaxially strained graphene. We observe the impact of the strain parameter in the emergence of towers of infinitely many bound states and the changes that it induces in the Efimov-like scaling of these cascades of bound states of the pristine case. We discuss our findings and present our conclusions in Section 4. An appendix with useful formulae is also included.

II. MONOPOLE IMPURITY AND A MAGNETIC FIELD IN STRAINED GRAPHENE
We start from the 2D Dirac Hamiltonian describing an electric monopole impurity in graphene immersed in a perpendicular magnetic field B uniformly distributed in space. Considering the tensor character of the Fermi velocity induced by strain, [35][36][37][38] , such an equation is written as where V (r ) corresponds to impurity potential under strain, a, b are the strain parameters along each spatial direction, π i = −i∂ i + (e/ c)A i with i = x, y, is the i-th component of the canonical momentum, A = B0 2 (−y, x) is the symmetric gauge vector potential describing the external magnetic field B, σ i are the Pauli matrices, and ∆ is the quasiparticle mass gap. The two component spinor Ψ ξs carries the valley (ξ = ±) and spin (s = ±) indices. Also, the standard convention Ψ +s = (ψ A , ψ B ) K T +s whereas Ψ −s = (ψ B , ψ A ) K T −s is used. Here, A, B refer to the two sublattices of a hexagonal graphene crystal. Because the interaction V (r ) does not depend on spin, this quantum number is omitted throughout. Then, by defining the parameters the equation governing the system can be written as where we have defined the operators and the quantities Here, l B = c/|eB 0 | is the magnetic length. Defining the following azimuthal-symmetric relations which describe the ellipse we have that hence obtaining the following form of the operators in (4), with the substitution ρ = r/l B in the last expression.
For further details about how the dimensionless quantities a, b (Eq. (2)) depend on the lattice parameters and the tensile strain ε, that measures compress and tensile deformations on the layer, the reader is referred to 35,38 .
In addition, the transformation showed in Eq. (6) allows us to solve the problem in a new coordinate system in which the ζ-dependency remains implicitly. More precisely, r 2 = ζ −1 x 2 + ζ y 2 and tan(φ) = ζ y/x. On the other hand, considering strain in the Coulomb potential (q/r ), with r 2 = ζr 2 cos 2 φ + ζ −1 r 2 sin 2 φ from Eq. (7), we found that thus, obtaining the zero-order approximation of the strained Coulomb potential with the charge redefinition q = qζ −1/2 . Unfortunately, it is well known that the problem of the Dirac equation in a magnetic field plus a Coulomb potential cannot be solved exactly 34 . Thus, considering a finite-size monopole impurity and neglecting the Coulombic tail (to avoid the difficulties of the analytical examination of this problem), we regularize the Coulomb potential at the origin by considering a potential well V (r ) = −V 0 θ(r 0 −r ), including the effects of strain in the radial variable r. The argument of the Heaviside function Using Eq. (56) of the Appendix, and noting that (13) so that, the potential well remains invariant under strain. Then, from Eq. (3) we get where (15) The spinor Ψ(r) can be written in terms of the eigenfunctions of the conserved angular momentum J z = L z + σ z /2 = −i∂ φ + σ z /2 as follows: with j = ±1/2, ±3/2,... Substituting this spinor in Eq. (14), through the replacement x = ρ 2 /2, we have for the upper component (17) which correspond to a Whittaker differential equation, whose general solution is expressed in terms of confluent hypergeometric functions M (a, c; x) and U (a, c; x) 39 as From this solution, we can obtain the general expression for g(ρ) by acting with S + upon the upper component of the spinor ψ A . Moreover, from the asymptotic behavior of M (a, c; x) and U (a, c; x), the regular solutions at the origin ρ → 0 (int) and ρ → ∞ (ext) for each component are, respectively, where C 1 and C 2 are arbitrary constants, Γ(n) denotes the Gamma function 39 and c = j + 1/2. These expressions are valid for j = ±1/2, ±3/2, ±5/2, . . . Continuity of the spinors at the boundary imply that condition that prescribes a transcendental relation for the energy eigenvalues of the system with total angular momentum j. Explicitly: This equation reduces to the corresponding relation for the unstrained case considered in 34 when λ = 1. In fact, a straightforward analysis follows merely replacing v F → λv F . Notice that because the functions M (a, c; x) and U (a, c; x) are not well defined for c = 0, −1, −2, . . ., which in our case correspond to negative values of j, by using the identities (58) and (59) of the Appendix, we are lead to the alternative expression valid for j ≤ −1/2: which we analyze below in some relevant regimes.

A. Collapse without external magnetic field
Let us first analyze the case of the potential well in graphene without the external magnetic field (B = 0 ⇒ l B → ∞ ⇒ p 2 v → ∞). By using the identities in Eq. (62) of the Appendix and taking z = p 2 v ρ 2 0 /2 ⇒ ρ 2 0 /2 = z/p 2 v , Eq. (24) becomes Without loss of generality, let us focus on valley K − (ξ = −1) 40 . The energy spectrum consists of a continuum part for |E| > ∆, and a discrete one for |E| < ∆, corresponding to bound states. The first bound state, E g ∆, is associated to the smallest centrifugal barrier j = −1/2, and appears at an arbitrarily small potential strength V 0 . Under these considerations, the quantities β and β in Eq. (27) where γ E is the Euler-Mascheroni constant (γ E ≈ 0.5772). The dependence of the ground state energy with respect to size of the monopole impurity r 0 at fixed V 0 is shown in Fig. 1. In the unstrained case, as the size of the impurity increases, the ground-state energy lowers till r 0 reaches a critical size where E g reaches E g = −∆, causing collapse. When we turn on the strain such that λ < 1, collapse occurs for impurities of smaller size. The quoted value λ =0.9 is consistent with the reported values of strain for which the graphene membrane can be contracted before folding 41 . According to 35,38 , this value for λ would obtain by tensile deformations around ε = 10% either along the zigzag direction or the armchair direction. For λ > 1, impurities need to be larger in order for strain to drive the collapse. The value λ =1.05 corresponds to compression deformation of 10% either along the zigzag direction or the armchair direction. It is worth to mention that the maximum reported value of strain before the crystal breaks corresponds to λ = 1.25 42 . On the other hand, as V 0 grows, the ground state gets closer to the threshold E = −∆. In Fig. 2, we show the surface of the parameters λv F , r 0 and V 0 /∆ for which E g = −∆. Crossing occurs at the critical value of V cr 0 , where j 0,1 ≈ 2.41 is the first zero of the function J 0 (x) 43 . From Eq. (29), it can be seen that values of λ > 1 lead to a higher value of V cr 0 , whereas for λ < 1, this threshold is smaller. We can think of this observation, as the strain deforms the lattice, the Dirac cones have now an elliptical cross section. Thus, the electrons orbiting near the larger semi-axis last more in their spiral-falling trajectory than those that are near the lower semi-axis.

B. Collapse including magnetic field
Let us now switch-on the effect of the external magnetic field. First, notice that neglecting the potential well potential, Eq. (24) in the K point and with j ≤ −1/2 render the energy eigenvalues in the well-known form of LLs corrected by strain through the renormalization of the Fermi velocity v F → λv F 35-38 , namely, and j + 1/2 ≤ n. It is important to mention that these states are degenerate because the energy only depends on the principal number n and not on the angular momentum value j. This degeneracy is lifted when a nonzero potential V 0 is considered. In Fig 3, we plot the energy eigenvalues for j = −1/2 and j = −3/2 for the extreme values of strain λ =0.9 and λ =1.05, both for the valley ξ = −. There are some important facts to notice: first, as we said before, the states with different angular momentum are no longer degenerate, and more and more of them decay from E = ∆ to E = −∆ as V 0 grows. This represents the fact that these states are diving into the (artificial) nucleus, taking place the atomic collapse. In addition, comparing with the case of unstrained graphene 34 , an expansive strain inhibits the atomic collapse from appearing, whereas an contracting one promotes collapse. The first crossing to the continuum of energy is given by the state with j = −1/2. The critical value of the potential V cr 0 in which this occur (when E g = −∆) can be obtained from Eq. (24), by noting that in this case Making use of U (0, 0; z) = U (0, 1; z) = 1, we have where a = l 2 B V cr 0 (V cr 0 −2∆)/2( λv F ) 2 . The critical value potential as a function of the quasiparticle mass gap is depicted in Fig. 4 for different values of the parameters ρ and λ. It is clearly seen that V cr 0 decreases drastically as the magnetic field strength B 0 increases (because ρ 0 = r 0 /l B and l B ∼ 1/|B 0 |) for fixed values of r 0 and ∆.
From the physical point of view, it represents the fact that the external magnetic field B 0 brings electronic orbitals closer to the nucleus, and thus, the attraction between an orbiting electron and a nucleus gets more intense, lowering the critical value of the potential. In addition, the value of the parameter λ, which contains information about the strain, provokes the critical po-tential value V cr 0 to grow for λ > 1 or to decrease for λ < 1 (Fig. 4), for fixed values of r 0 , ∆ and B 0 . This result leads us to conclude the possibility to induce such an strain in the graphene lattice which makes electronic collapse easier to take place in one direction and more difficult in the other, being capable to control an effective directed-current. Another observation is that in all the cases with B 0 = 0, V 0,cr goes to zero as ∆ decreases arbitrarily. From that, we can conclude that the presence of an homogeneous magnetic field catalyzes atomic collapse in the Coulomb center doped graphene for any strain. In particular, in the case of gapless graphene in the presence of an homogeneous magnetic field atomic collapse would take place for any value of Z.

III. POINT DIPOLE IN STRAINED GRAPHENE
In this section we address the case of a point electric dipole placed in a graphene plane 19 , but considering additionally the effect of uniform uniaxial strain. The electric dipole consists in two opposite charges +Q and −Q separated a distance d. Thus, its potential taking strain into account is 44 where p = Qd is the dipole moment, and are the lengths of the vectors r + , r − under strain joining the negative and positive charge positions with the observation point respectively. From Eq. (6) with c 1 = ζ −2 − 1 and c 2 = d/2. So that and thus, in the limit d → 0, the point electric dipole potential under strain is with the redefinition of the dipole moment p = ζ −1 p.
The corresponding Hamiltonian has the form where we are introducing the strain effects by the tensor Fermi velocity as in the previous section. This Hamiltonian has an intrinsic electron-hole symmetry U HU † = −H, with the unitary operator U = σ x R x satisfying U 2 = 1, where R x is the operator of reflection x → −x. An eigenstate Ψ E (x, y) with energy E is mapped to another eigenstate with energy −E as Hence, all solutions of the Dirac equation come in pairs with ±E, and it is enough to study one of the eigenstates to automatically know the other. Thus, the Dirac equation governing the system is written as (39) where P x = −i∂ x and P y = −i∂ y are the x, y components of the linear momentum, and Ψ = (ψ A ψ B ) T . Thus, the Dirac equation turns into and thus we are lead to the equation for ψ B of the form where ∇ 2 = ∂ 2 r + 1 r ∂ r + 1 r 2 ∂ 2 φ is the 2D Laplace operator in cylindrical coordinates. Let us propose a separable solution of the form which upon separation of variables can be written as where γ is the separation constant. The angular equation which is independent of the parameter , corresponds to a Mathieu equation 43,45 , whose solutions where κ = ± is the parity, are 2π-periodic and are only permitted for characteristic values of γ, where the angular momentum j = 0, 1, 2, . . . is unconventional because of the anisotropy, satisfying the relation j + κ ≥ 0. These characteristic values are ordered for a given value of p as γ 0,+ (p ) < γ 1,− (p ) < γ 1,+ (p ) < γ 2,− (p ) < ... The radial equation is given by and has the form of a MacDonald equation 45 whose solutions are the modified Bessel functions K ν (z), which in our case correspond to ν = √ γ and z = (r √ 2∆ )/( λv F ). Recall that for bound states ( > 0), these solutions have to decay for r → 0 and to be regular at the origin. For this purpose, we use the regularization condition R(r 0 ) = 0, which gives the energy quantization condition within each tower of (j, κ), where z n are the zeros of the function K √ γj,κ (z) satisfying z 1 < z 2 < z 3 < ... Because this function only has zeros for √ γ j,κ imaginary 45 , then the condition γ j,κ (p ) < 0 is required for bound states. This is satisfied for values p > p j,κ with γ j,κ (p j,κ ) = 0.
These values of p j,κ have been already presented in 19 .
By increasing the value of p , each time that p hits a critical value p j,κ , a new infinite tower of bound states emerges from the continuum. Then, by expanding K is (z) for small z as in Eq. (64) (with s j,κ (p ) = −γ j,κ (p ) for p > p j,κ ), we obtain an explicit expression for the bound energies near the edges ±∆ given by n,j,κ = where φ(s) = (2/s) arg Γ(1+is), and the values s j,κ are 19 where α ≈ 0.956. It can be seen that these states verify an Efimov-like universal scaling law as in 19 n+1 n = e −2π/sj,κ , but with the λ-dependence inserted in the s j,κ terms. Plots of the energy spectrum versus strained dipole moment p for different values of the strain strength λ are presented in Fig. 5. In all of them the first electron, n,j,κ , bound states are presented for each tower, decaying from the continuum to discrete values. Just as in the case of a single charged impurity seen in the previous section, the role of the strain strength is to promote (for values of λ < 1) or inhibit (for values of λ > 1) the decay of the electron and hole bound states, making the collapse in both systems easier or more difficult to occur depending on the strain. Also, for n → ∞ we have that → 0, so that the edges ±∆ are accumulation points. These facts suggest that electrons can be captured by a dipole potential in graphene, and strain can be used to tune the strength of confinement and to promote or inhibit atomic collapse.

IV. CONCLUSIONS
In this article we have carried out an analysis of the bound states for the anisotropic 2D Dirac equation (by virtue of uniform uniaxial strain) in the potential of point electric monopole and dipole impurities in the context of graphene physics. In the former case, in absence of strain, it is known that the combined effect of a uniform magnetic field and a Coulomb impurity yield to a scenario of atomic collapse. Analytically, this problem cannot been solved. Nevertheless, by regularizing the Coulomb potential at the origin through a potential well, the problem exhibits interesting features, like atomic collapse in absence of magnetic fields provided the potential well has a value above certain critical number strongly connected with the size of the impurity. The magnetic field drives this critical strength to zero, such that atomic collapse occurs for impurities of arbitrary size. Furthermore, as the magnetic field increases its strength, more and more states with angular momentum j < −1/2 start to dive into the continuum. On the other hand, for the point dipole impurity, even though no collapse of states is observed, the appearance of cascades of infinitely many bound states with a universal Efimov-like scaling appears as a reminiscent of the true collapse that happens in the case of a finite-size dipole. Because the role of strain is seen in a renormalization of the Fermi velocity v F → λv F , for λ < 1 these two scenarios are promoted to occur under less restrictive conditions, as compared with the ideal case. For λ > 1, the situation is the opposite. As an additional result, we observe that the effect of strain in the Coulomb and point electric dipole potentials addressed in this work can be understood as a redefinition of the electric charge and the dipole moment, respectively, in terms of the strain parameter ζ.
The role of position-dependent strain is currently under consideration in our group, with the addition of the generation of electric fields due to this deformations. Along these lines, we point out the recent findings regarding the collapse of Landau levels under strain and a uniform electric field recently discussed in 46 . Results will be presented elsewhere.