Bound States and Supercriticality in Graphene-Based Topological Insulators

We study the bound state spectrum and the conditions for entering a supercritical regime in graphene with strong intrinsic and Rashba spin-orbit interactions within the topological insulator phase. Explicit results are provided for a disk-shaped potential well and for the Coulomb center problem


Introduction
The electronic properties of graphene monolayers are presently under intense study. Previous works have already revealed many novel and fundamental insights; for reviews, see [1,2]. Following the seminal work of Kane and Mele [3], it may be possible to engineer a two-dimensional (2D) topological insulator (TI) phase [4] in graphene by enhancing the-usually very weak [5][6][7]-spin-orbit interaction (SOI) in graphene. This enhancement could, for instance, be achieved by the deposition of suitable adatoms [8].
Remarkably, random deposition should already be sufficient to reach the TI phase [9][10][11] where the effective "intrinsic" SOI ∆ exceeds (half of) the "Rashba" SOI λ. So far, the only 2D TIs realized experimentally are based on the mercury telluride class. Using graphene as a TI material constitutes a very attractive option because of the ready availability of high-quality graphene samples [1] and the exciting prospects for stable and robust TI-based devices [4], see also [12,13].
In this work, we study bound-state solutions and the conditions for supercriticality in a graphene-based TI. Such questions can arise in the presence of an electrostatically generated potential well ("quantum dot") or for a Coulomb center. The latter case can be realized by artificial alignment of Co trimers [14], or when defects or charged impurities reside in the graphene layer. Without SOI, the Coulomb impurity problem in graphene has been theoretically studied in depth [15][16][17][18][19][20]; for reviews, see [1,2]. Moreover, for λ = 0, an additional mass term in the Hamiltonian corresponds to the intrinsic SOI ∆ (see below), and the massive Coulomb impurity problem in graphene has been analyzed in [21][22][23][24][25][26]. However, a finite Rashba SOI λ is inevitable in practice and has profound consequences. In particular, λ = 0 breaks electron-hole symmetry and modifies the structure of the vacuum. We therefore address the general case with both ∆ and λ finite, but within the TI phase ∆ > λ/2, in this paper. Experimental progress on the observation of Dirac quasiparticles near a Coulomb impurity in graphene was also reported very recently [14], and we are confident that the topological version with enhanced SOI can be studied experimentally in the near future. Our work may also be helpful in the understanding of spin-orbit mediated spin relaxation in graphene [27].
The atomic collapse problem for Dirac fermions in an attractive Coulomb potential, V (r) = −hv F α/r, could thereby be realized in topological graphene. Here we use the dimensionless impurity strength where Z is the number of positive charges held by the impurity; κ a dielectric constant characterizing the environment; and v F ≈ 10 6 m/s the Fermi velocity. Without SOI, the Hamiltonian is not self-adjoint for α > α c = 1/2, and the potential needs short-distance regularization, e.g., by setting V (r < R) = −hv F α/R with short-distance cutoff R of the order of the lattice constant of graphene [1,2]. Including a finite "mass" ∆, i.e., the intrinsic SOI, but keeping λ = 0, the critical coupling α c is shifted to [24] α c 1 2 + π 2 ln 2 (0.21∆R/hv F ) approaching the value α c = 1/2 for R → 0. In the supercritical regime α > α c , the lowest bound state "dives" into the valence band continuum (Dirac sea). It then becomes a resonance with complex energy, where the imaginary part corresponds to the finite decay rate into the continuum. Below we show that the Rashba SOI provides an interesting twist to this supercriticality story. The structure of this article is as follows. In Section 2 we introduce the model and summarize its symmetries. The case of a circular potential well is addressed in Section 3 before turning to the Coulomb center in Section 4. Some conclusions are offered in Section 5. Note that we do not include a magnetic field (see, e.g., [28,29]) and thus our model enjoys time-reversal symmetry. Below, we often use units withh = v F = 1.

Kane-Mele Model with Radially Symmetric Potential
We study the Kane-Mele model for a 2D graphene monolayer with both intrinsic (∆) and Rashba (λ) SOI [3] in the presence of a radially symmetric scalar potential V (r). Assuming that V (r) is sufficiently smooth to allow for the neglect of inter-valley scattering, the low-energy Hamiltonian near the K point (τ = +1) is given by with Pauli matrices σ x,y,z (s x,y,z ) in sublattice (spin) space [1]. The Hamiltonian near the other valley (K point) follows for τ = −1 in Equation (3). We note that a sign change of the Rashba SOI, λ → −λ, does not affect the spectrum due to the relation H(−λ) = s z H(λ)s z . Without loss of generality, we then put ∆ ≥ 0 and λ ≥ 0.
Using polar coordinates, it is now straightforward to verify (see also [21]) that total angular momentum, defined as is conserved and has integer eigenvalues j. For given j, eigenfunctions of H must then be of the form Next we combine the radial functions to (normalized) four-spinors In this representation, the radial Dirac equation for total angular momentum j and valley index with Hermitian matrix operators (note that ∆ denotes the intrinsic SOI and not the Laplacian) where we use the notation ∇ One easily checks that Equation (8) satisfies the parity symmetry relation Note that this "parity" operation for the radial Hamiltonian is non-standard in the sense that the valley is not changed by the transformation σ y s y , spin and sublattice are flipped simultaneously, and only the y-coordinate is reversed. (We will nonetheless refer to σ y s y as parity transformation below.) A second symmetry relation connects both valleys, Using Equation (10), this relation can be traced back to a time-reversal operation. Equations (10) and (11) suggest that eigenenergies typically are four-fold degenerate.
When projected to the subspace of fixed (integer) total angular momentum j, the current density operator has angular component J φ = σ x and radial component J r = −τ σ y for arbitrary j. When real-valued entries can be chosen in Φ j,τ (r), the radial current density thus vanishes separately in each valley. We define the (angular) spin current density as J S φ = s z σ x . Remarkably, the transformation defined in Equation (11) conserves both (total and spin) angular currents, while the transformation in Equation (10) reverses the total current but conserves the spin current. Therefore, at any energy, eigenstates supporting spin-filtered counterpropagating currents are possible. However, in contrast to the edge states found in a ribbon geometry [3], these spin-filtered states do not necessarily have a topological origin.
We focus on one K point (τ = +) and omit the τ -index henceforth; the degenerate τ = − Kramers partner easily follows using Equation (11). In addition, using the symmetry (10), it is sufficient to study the model for fixed total angular momentum j ≥ 0.

Zero Total Angular Momentum
For arbitrary V (r), we now show that a drastic simplification is possible for total angular momentum j = 0, which can even allow for an exact solution. Although the lowest-lying bound states for the potentials in Sections 3 and 4 are found in the j = 1 sector, exact statements about what happens for j = 0 are valuable and can be explored along the route sketched here.
The reason why j = 0 is special can be seen from the parity symmetry relation in Equation (10). The parity transformation σ y s y connects the ±j sectors, but represents a discrete symmetry of the j = 0 radial Hamiltonian H j=0,τ [see Equation (8)] acting on the four-spinors in Equation (6). Therefore, the j = 0 subspace can be decomposed into two orthogonal subspaces corresponding to the two distinct eigenvalues of the Hermitian operator σ y s y . This operator is diagonalized by the matrix In fact, using this transformation matrix to carry out a similarity transformation,H j, For j = 0, the upper and lower 2 × 2 blocks decouple. Each block has the signature ("parity") σ = ± corresponding to the eigenvalues in Equation (13), and represents a mixed sublattice-spin state, see Equations (6) and (12). For parity σ = ±, the 2 × 2 block matrix in Equation (14) is formally identical to an effective λ = 0 problem with j = 0, fixed s z = σ, and the substitutions This implies that for j = 0 and arbitrary V (r), the complete spectral information for the full Kane-Mele problem (with λ = 0) directly follows from the λ = 0 solution.

Solution in Region with Constant Potential
We start our analysis of the Hamiltonian (3) with the general solution of Equation (7) for a region of constant potential. Here, it suffices to study V (r) = 0, since E and V enter only through the combination E − V in Equation (8). In Section 3, we will use this solution to solve the case of a step potential.
The general solution to Equation (7) follows from the Ansatz where the c i are real coefficients; B j is one of the cylinder (Bessel) functions; j ; and p denotes a real spectral parameter. In particular, √ p is a generalized radial wavenumber. We here assume true bound-state solutions with real-valued energy. However, for quasi-stationary resonance states with complex energy, p and the c i may be complex as well.
Using the Bessel function recurrence relation, ∇ , the set of four coupled differential Equations (7) simplifies to a set of algebraic equations Notably, j does not appear here, and therefore the spectral parameter p depends only on the energy E. The condition of vanishing determinant then yields a quadratic equation for p, with the two solutions Which Bessel function is chosen in Equation (16) now depends on the sign of p ± and on the imposed regularity conditions for r → 0 and/or r → ∞.
For p ± > 0, a solution regular at the origin is obtained by putting B j = J j , which describes standing radial waves. Equation (17) then yields the unnormalized spinor For p ± < 0, instead it is convenient to set B j = H (1) j in Equation (16). Using the identity H , the unnormalized spinor resulting from Equation (17) then takes the form where the modified Bessel function K j ( √ −p ± r) describes evanescent modes, exponentially decaying at infinity.

Solution without Potential
In a free system, i.e., when V (r) = 0 for all r, the only acceptable solution corresponding to a physical state is obtained when p ± > 0 [30]. For ∆ < λ/2, at least one p ± > 0 in Equation (18) for all E, and the system is gapless. However, the TI phase defined by ∆ > λ/2 has a gap as we show now.
For ∆ > λ/2, Equation (18) tells us that for E > ∆ and for E < E − , both solutions p ± are positive and hence (for given j and τ ) there are two eigenstates Φ j,p ± for given energy E. However, within the energy window [with E ± in Equation (18)] we have p + > 0 and p − < 0, i.e., only the eigenstate Φ j,p + represents a physical solution. Both p ± are negative when E + < E < ∆, and no physical state exists at all. This precisely corresponds to the topological gap in the TI phase [3]. Note that due to the Rashba SOI, the valence band edge is characterized by the two energies E ± , with halved density of states in the energy window (21). One may then ask at which energy (E + or E − ) the supercritical diving of a bound state level in an impurity potential takes place.

Bound States
In this section, we study a circular potential well with radius R and depth V 0 > 0 We always stay within the TI phase ∆ > λ/2, where bound states are expected for energies E = E B in the window max(∆ − V 0 , E + ) < E B < ∆. For r < R, the corresponding radial eigenspinor [see Equation (6)] is written with arbitrary prefactors A < ± in the form with Equation (19) for Φ j,p ± (r). Here, thep ± > 0 follow from Equation (18) by including the potential shift,p For r > R, the general solution is again written as However, now Φ j,p ± is given by Equation (20), since p ± < 0 for true bound states with only evanescent states outside the potential well. The continuity condition for the four-spinor at the potential step, Φ < j (R) = Φ > j (R), then yields a homogeneous linear system of equations for the four parameters (A <,> ± ). A nontrivial solution is only possible when the determinant of the corresponding 4 × 4 matrix C(E) (which is too lengthy to be given here but follows directly from the above expressions) vanishes Solving the energy quantization condition (26) then yields the discrete bound-state spectrum (E B ). It is then straightforward to determine the corresponding spinor wavefunctions. Numerical solution of Equation (26) yields the bound-state spectrum shown in Figure 1. When V 0 exceeds a (j-dependent) "threshold" value, V t , a bound state splits off the conduction band edge. When increasing V 0 further, this bound-state energy level moves down almost linearly, cf. inset of Figure 1, and finally reaches the valence band edge E + = −∆ + λ at some "critical" value V 0 = V c . (For j = 0, we will see below that this definition needs some revision.) Increasing V 0 even further, the bound state is then expected to dive into the valence band and become a finite-width supercritical resonance, i.e., the energy would then acquire an imaginary part.

Zero Angular Momentum States
Surprisingly, for j = 0, we find a different scenario where supercritical diving, with finite lifetime of the resonance, happens only for half of the bound states entering the energy window (21). Noting that states with different parity σ = ± do not mix, see Section 2.2, we observe that all σ = + bound states enter the valence band as true bound states (no imaginary part) throughout the energy window (21) while the valence band continuum is spanned by the σ = − states. We then define V c for (j = 0, σ = +) bound states as the true supercritical threshold where E B = E − = −∆ − λ. However, the (j = 0, σ = −) bound states become supercritical already when reaching E + = −∆ + λ.
Therefore an intriguing physical situation arises for j = 0 in the energy window (21). While σ = + states are true bound states (no lifetime broadening), they coexist with σ = − states which span the valence band continuum or possibly form supercritical resonances. For E < E − , however, all bound states dive, become finite-width resonances, and eventually become dissolved in the continuum.

Threshold for Bound States
Returning to arbitrary total angular momentum j, we observe that whenever V 0 hits a possible threshold value V t , a new bound state is generated, which then dives into the valence band at another potential depth V 0 = V c (and so on). Analytical results for all possible threshold values V t follow by expanding Equation (26) for weak dimensionless binding energy δ ≡ 1 − E B /∆. For δ 1 and j = 1, Equation (26) yields after some algebra where γ ≈ 0.577 is the Euler constant andλ = λ/2∆. The binding energy approaches zero for V 0 → 0, where Equation (27) simplifies to For vanishing Rashba SOI λ = 0, this reproduces known results [25]. For any λ < 2∆, we observe that the j = 1 bound state in Equation (28) exists for arbitrarily shallow potential depth V 0 .
The threshold values V t for higher-lying j = 1 bound states also follow from the binding energy (27), since δ vanishes for J 1 (z + ) = 0 and for J 1 (z − ) = 0. When one of these two conditions is fulfilled at some V 0 = V t , a new bound state appears for potential depth above V t . This statement is in fact quite general: By similar reasoning, we find that the threshold values V t for j = 0 follow by counting the zeroes of J 0 (z ± ). Without SOI, this has also been discussed in [31]. Note that this argument immediately implies that no bound state with j = 0 exists for V 0 → 0.
From the above equations, we can then infer the threshold values V t for all bound states with j = 0 or j = 1 in analytical form. These are labeled by n = 1, 2, . . . and σ = ± (for j = 0, σ corresponds to parity) V t,j,n,± = (∆ ± λ/2) −1 + 1 + γ 2 j,n /[R(∆ ± λ/2)] 2 (29) where γ j,n is the nth zero of the J j Bessel function. Likewise, for j > 1, the condition for the appearance of a new bound state is Close examination of this condition shows that no bound states with j > 1 exist for V 0 → 0. We conclude that bound states in a very weak potential well exist only for j = 1.

Supercritical Behavior
As can be seen in Figure 1, the lowest j = 1 bound state is also the first to enter the valence band continuum for V 0 = V c . For λ = 0, the critical value is known to be [25] V c = ∆ 1 + 1 + with γ 0,1 ≈ 2.41. The energy of the resonant state acquires an imaginary part for V 0 > V c [25]. For λ > 0, we have obtained implicit expressions for V c , plotted in Figure 2. Note that these results reproduce Equation (31) for λ → 0. The almost linear decrease of V c with increasing λ, see Figure 2, can be rationalized by noting that the valence band edge is located at E + = −∆ + λ. Thereby supercritical resonances could be reached already for lower potential depth by increasing the Rashba SOI. Similarly, with increasing disk radius R, the critical value V c decreases, see the inset of Figure 2. For the lowest (j = 0, σ = ±) bound state, the critical value in fact follows in analytical form, where γ 1,1 ≈ 3.83. Since the parity decoupling in Section 2.2 only holds for j = 0, it is natural to expect that all j = 0 bound states turn into finite-width resonances when E B < E + . This expectation is confirmed by an explicit calculation as follows. Within in the window E − < E B < E + , a true bound state should not receive a contribution from Φ j,p + >0 (r) for r > R, but instead has to be obtained by matching an Ansatz as in Equation (23) for the spinor state inside the disk (r < R) to an evanescent spinor state ∝ Φ j,p − <0 (r > R). However, the matching condition is then found to have no real solution E B , i.e., there are no true bound states with j = 0 in the energy window (21). We therefore conclude that all j = 0 bound states turn supercritical when E B < E + . Note that this statement includes the lowest-lying bound state (which has j = 1). This implies that a finite Rashba SOI can considerably lower the potential depth V c required for entering the supercritical regime.

Coulomb Center
We now turn to the Coulomb potential, V (r) = −α/r, generated by a positively charged impurity located at the origin, with the dimensionless coupling strength α in Equation (1). We consider only the TI phase ∆ > λ/2 and analyze the bound-state spectrum and conditions for supercriticality. Again, without loss of generality, we focus on the K point only (τ = +), and first summarize the known solution for λ = 0 [2,21,24]. In that case, s z = ± is conserved, and the spin-degenerate bound-state energies are labeled by the integer angular momentum j and a radial quantum number n = 1, 2, 3, . . . (for j > 0, n = 0 is also possible) The corresponding eigenstates then follow in terms of hypergeometric functions. The lowest bound state is E j=1,n=0 = ∆ √ 1 − 4α 2 , which dives when α = α c = 1/2; note that α c precisely corresponds to V c in Section 3. In particular, for (j = 0, σ = ±) states we define α c in the same manner. Next we discuss how this picture is modified when the Rashba coupling λ is included.
Following the arguments in Section 2.2 for j = 0, the combination of Equation (33) with Equation (15) immediately yields the exact bound-state energy spectrum (n = 1, 2, 3, . . .) E j=0,n,σ=± = (∆ ± λ/2) The corresponding eigenstates then also follow from [21,24]. The very same reasoning also applies to a regularized 1/r potential [23,24], where V (r < R) is replaced by the constant value V = −α/R. Here, R is a short-distance cutoff scale of the order of the lattice spacing. The solution of the bound-state problem then requires a wavefunction matching procedure, which has been carried out in [24]. Thereby we can already infer all bound states for j = 0. Figure 3 shows the resulting j = 0 bound-state spectrum vs. α for the regularized Coulomb potential. Within the energy window Equation (21), we again find that states with parity σ = + remain true bound states that dive only for E B < E − , while σ = − states show supercritical diving already for E B < E + . Figure 4 shows the corresponding critical couplings α c for σ = ±, where the lowest j = 0 bound state with parity σ turns supercritical. Note that for finite R and λ → 0, a unique value for α c is found, while for λ = 0 two different critical values for α c are found. However, this conclusion holds only for finite regularization parameter R, i.e., it is non-universal. As seen in the inset of Figure 4, in the limit R → 0, both critical values for α c approach α c = 1/2 again, which is the value found without SOI.
Finally, for j = 0, we can then draw the same qualitative conclusions as in Section 3.4 for the potential well. In particular, we expect that all j = 0 bound states turn supercritical when their energy E B reaches the continuum threshold at E B = E + = −∆ + λ.

Conclusions
In this work, we have analyzed the bound-state problem for the Kane-Mele model of graphene with intrinsic (∆) and Rashba (λ) spin-orbit couplings when a radially symmetric attractive potential V (r) is present. We have focused on the most interesting "topological insulator" phase with ∆ > λ/2. The Rashba term λ leads to a restructuring of the valence band, with a halving of the density of states in the window E − < E < E + , where E ± = −∆ ± λ. This has spectacular consequences for total angular momentum j = 0, where the problem can be decomposed into two independent parity sectors (σ = ±). The σ = + states remain true bound states even inside the above window and coexist with the continuum solutions as well as possible supercritical resonances in the σ = − sector. However, all j = 0 bound states exhibit supercritical diving for E < E + , where the critical threshold (V c or α c for the disk or the Coulomb problem, respectively) is lowered when the Rashba term is present. We hope that these results will soon be put to an experimental test.