Two-electron bound states near a Coulomb impurity in gapped graphene

# Two-electron bound states near a Coulomb impurity in gapped graphene

Alessandro De Martino Department of Mathematics, City, University of London, EC1V 0HB London, United Kingdom    Reinhold Egger Institut für Theoretische Physik, Heinrich-Heine-Universität, D-40225 Düsseldorf, Germany
August 4, 2019
###### Abstract

We formulate and solve the perhaps simplest two-body bound state problem for interacting Dirac fermions in two spatial dimensions. A two-body bound state is predicted for gapped graphene monolayers in the presence of weakly repulsive electron-electron interactions and a Coulomb impurity with charge , where the most interesting case corresponds to . We introduce a variational Chandrasekhar-Dirac spinor wave function and show the existence of at least one bound state. This state leaves clear signatures accessible by scanning tunneling microscopy. One may thereby obtain direct information about the strength of electron-electron interactions in graphene.

## I Introduction

Ever since the Nobel prize winning experiments by the Manchester group in 2004 geim2004 (); geim2005 (), two-dimensional (2D) graphene monolayers continue to attract a lot of attention. Besides graphene’s application potential, much of the interest comes from the fact that low-energy fermionic quasiparticles in graphene are governed by the 2D Dirac Hamiltonian rmp1 (); goerbig (); philip (); eva (); kotov (); miransky (). As a consequence, typical effects predicted by relativistic quantum mechanics and/or quantum electrodynamics can be studied in table-top experiments. Recent progress has demonstrated that one can reach the ballistic (disorder-free) regime Kim2016 (), e.g., by encapsulating the graphene layer in boron nitride crystals boron (). We will therefore not consider disorder effects in this work.

Here we address the perhaps simplest interacting problem for relativistic fermions by considering a gapped graphene monolayer subject to relatively weak (screened) Coulomb interactions and in the presence of a single Coulomb impurity. A band gap in the Dirac fermion spectrum may be caused by a variety of mechanisms, e.g., by strain engineering vozmediano (), spin-orbit coupling kanemele2005 (); huertas (), substrate-induced superlattices ponomarenko (); song (), or simply as finite-size effect in a ribbon geometry rmp1 (). The Coulomb interaction strength is customarily quantified in terms of the effective fine-structure constant

 α=e2κℏvF≃2.2κ, (1)

where is a dielectric constant due to the surrounding substrate and ms denotes the Fermi velocity kotov (). The estimate in Eq. (1) follows with and . The value of thus depends on the choice of gating geometry and the ensuing screening effects. As explained below, the weak interaction regime where our theory is applicable is defined by . In this regime, we find that a pair of repulsively interacting Dirac fermions subject to the attractive potential of a Coulomb impurity with charge will form a two-body bound state localized near the impurity. We mainly focus on the most interesting case , but also comment on what happens for .

The signatures of the predicted bound state could be observed by means of scanning tunnelling microscopy (STM) experiments similar to those previously reporting supercritical behavior in graphene crommie1 (); luican (); crommie2 () and trapped electron states in electrostatically defined graphene dots Lee2016 (); Gutierrez (). An experimental observation of the predicted two-body bound state could therefore probe and quantify Coulomb interaction effects. In order to consistently formulate this relativistic quantum-mechanical two-particle problem, it is necessary to stay away from the supercritical instability rmp1 (); kotov (); novikov (); gamayun (), otherwise one inevitably has to face a complicated many-body problem kotov (). For small values of , supercriticality is absent, and as explained below, a two-particle bound state exists. We mention in passing that the two-body problem in graphene has also been studied in Refs. guinea (); terekhov (); koch (); ziegler (); entin (); gaul (); shytov (); portnoi (). In contrast to our work, however, those papers considered translationally invariant settings without Coulomb impurity. Similar problems (again without Coulomb impurity) have also been analyzed in the context of Dirac excitons in single-layer transition-metal dichalcogenides rodin (); zhong (); trushin ().

The negatively charged two-electron hydrogen ion represents a classic problem of nonrelativistic quantum mechanics chandrasekhar (); bethe (); hill (); bransden (); andersen (); sorba (). In particular, it is well known hill () that has a single bound state in three spatial dimensions. The simplest way to prove its stability is to demonstrate the existence of a variational wave function with energy below the ground-state energy of the hydrogen atom. Interestingly, it is impossible to achieve this task with a factorized wave function sorba (). The simplest way to construct a variational wave function for the ground state of is due to Chandrasekhar chandrasekhar (). With denoting the distance of the respective electron from the proton, the Chandrasekhar ansatz for the spatial part of the two-particle wave function, , contains two variational parameters () and is (up to a normalization factor) given by

 Ψ(r1,r2)=e−ar1−br2−ϵe−br1−ar2, (2)

where corresponds to a spin singlet/triplet state, respectively. Here the important insight is that and are not required to be identical. Indeed, the minimal variational energy for a two-body bound state is obtained for in the spin singlet configuration (). Improved variational energy estimates can be obtained by taking into account the dependence on the relative distance in Eq. (2), e.g., through a correlation factor of the form sorba (). Since the inclusion of such a factor is technically cumbersome yet not expected to cause qualitative changes, cf. Ref. sorba (), we here restrict ourselves to uncorrelated wave functions and leave the analysis of correlation effects to future work.

The nonrelativistic hydrogen ion has also been studied for the 2D case. For instance, the so-called problem describes a donor impurity ion with two attached electrons in a 2D semiconductor quantum well phelps (); louie (); larsen92a (); larsen92b (); sandler (); ivanov (). Interesting experimental results on the problem have appeared in Ref. etienne (), where effects of quantum confinement on two-body bound-state energies have been studied. In the absence of a magnetic field, only a single bound state in the spin singlet sector is expected again. We note in passing that the problem is also similar to the negatively charged exciton () problem, which was experimentally studied in quantum wells pepper ().

In this paper, we turn to the 2D relativistic counterpart of the above system, which is realizable in gapped graphene monolayers (or topological insulator surfaces) containing a Coulomb impurity. The corresponding relativistic problem is more difficult to define (and solve) because the single-particle Dirac Hamiltonian is unbounded from below brown (); talman (); nakatsuji (). It then appears at first sight as if two-particle states of arbitrarily low energy can be generated by electron-electron interactions. As discussed below, in order to avoid this spurious and ultimately unphysical effect, it is necessary to employ a projection scheme which defines a mathematically clean framework. Such a projection scheme can be devised for interacting Dirac fermions in graphene if (i) a single-particle gap exists (), and (ii) electron-electron interactions are weak, see Refs. sucher1 (); sucher2 (); haeusleregger (); Egger ().

The structure of the remainder of this article is as follows. We introduce the Dirac-Coulomb model and the appropriate projection scheme in Sec. II, followed by a discussion of the variational approach to the two-body bound state problem in Sec. III. (Details have been delegated to two appendices.) To that end, we formulate a Chandrasekhar-Dirac spinor ansatz generalizing Eq. (2) to the relativistic case. We evaluate all needed matrix elements and discuss the variational bound-state energy as a function of . In particular, we show variational estimates for the energy of the bound state in the presence of a standard Dirac mass term causing a band gap, assuming a spin-singlet state and impurity charge . We then turn to generalizations in Sec. IV, where we shall address (i) what happens for impurity charge , (ii) the effects of a topologically nontrivial band gap as obtained, e.g., from spin-orbit coupling effects, and (iii) the role of the valley state of each quasiparticle. In Sec. V, we address the observable consequences of the two-body bound state accessible to STM experiments. Finally, we offer some conclusions in Sec. VI.

## Ii Dirac-Coulomb Model

We model the interacting two-particle problem for a gapped graphene monolayer in the presence of a charge impurity by the (properly projected, see below) Dirac-Coulomb Hamiltonian

 H=∑l=1,2HD(l)+V2b. (3)

Here is the usual single-particle massive Dirac-Weyl Hamiltonian for particle ,

 HD(l)=H0(l)+Hgap(l)+V1b(l), (4)

with kinetic part (the index being understood) rmp1 ()

 H0=vF(τzσxpx+σypy), (5)

where Pauli matrices () act in valley (sublattice) space and the momentum operator is . With the dielectric constant , see Eq. (1), the single-particle potential

 V1b=−Ze2κr=−ZαℏvFr (6)

describes a charge- impurity at the origin. We mainly focus on the most interesting case of unit charge, , but comment on the case in Sec. IV.1. In , we also include the term

 Hgap=Δσz, (7)

which results in a topologically trivial band gap of size rmp1 (). However, it is straightforward to generalize our analysis to the case of a spin- and valley-dependent topological band gap, e.g., due to an intrinsic spin-orbit coupling mechanism kanemele2005 (); huertas (),

 Hso=Δσzτzsz, (8)

where are Pauli matrices in spin space, see Sec. IV.2. Note that has the same sign for both spin projections while has opposite sign.

The operator is Hermitian only for . Indeed, the lowest single-particle bound-state energy becomes imaginary for , see Eq. (57) in App. A. We note that by regularizing the singularity of the Coulomb potential, the threshold shifts to a larger value , which weakly depends on the precise regularization prescription gamayun (). Even in the regularized scheme, however, Hermiticity is lost for . In the latter regime, bound states of the regularized Hamiltonian dive into the negative part of the continuum spectrum () one by one, and thereby become quasi-stationary states, see Ref. gamayun (). As remarked in Sec. I, one then necessarily has to consider the full many-body problem for . For this reason, we will restrict ourselves to the weak-coupling regime. Moreover, for the sake of transparency, we focus on the narrower range , where no need for regularization arises and the exact Dirac-Coulomb wave functions summarized in App. A can be used. For the most interesting case , this implies that our theory is at best applicable for .

Since is diagonal in valley and spin space, we can always choose a factorized form for the single-particle wave functions,

 Ψτ,s(x,y)=Ψ(x,y)|τ,s⟩, (9)

where refers to the spatial part (including sublattice space), and with and denotes the eigenstates of the operator . For an eigenstate of with energy , the symmetry relations

 σyΨ∗τ,s=Ψτ,−s,σyΨτ,s=Ψ−τ,s,Ψ∗τ,s=Ψ−τ,−s, (10)

are a manifestation of the well-known fourfold spin-valley degeneracy of each energy level rmp1 ().

Next, the two-body Coulomb interaction in Eq. (3) is formally given by (see App. B)

 V2b=e2κ|r1−r2|. (11)

However, since the spectrum of the single-particle Dirac Hamiltonian is unbounded from below, the many-body Dirac-Coulomb Hamiltonian (3) cannot have true two-particle bound states. This problem was pointed out long ago by Brown and Ravenhall in their study of relativistic effects in the helium atom brown ().

A rigorous procedure to construct a mathematically well-defined many-particle Hamiltonian for interacting massive Dirac fermions in three spatial dimensions was devised by Sucher sucher1 (); sucher2 (). Fortunately, as has been shown in detail in Refs. haeusleregger (); Egger (), Sucher’s approach can be readily adapted to the case of 2D Dirac fermions in graphene as long as the single-particle spectrum exhibits a band gap. In effect, in Eq. (3) thereby has to replaced by the projected Hamiltonian haeusleregger (); Egger ()

 H+=HD(1)+HD(2)+Λ+V2bΛ+, (12)

with the projection operator . Here, using , the single-particle operator projects onto the space spanned by the positive-energy eigenstates of . As detailed in Refs. sucher1 (); sucher2 (); haeusleregger (); Egger (), the projected Hamiltonian takes into account the most important effects of the electron-electron interaction. In fact, due to the presence of a band gap, the replacement does not introduce approximations concerning the ground state of the system in the limit of weak Coulomb repulsion. Moreover, the projection guarantees that the Hamiltonian can possess bona-fide two-particle bound states. For instance, in the nonrelativistic limit realized for energies very close to the upper band edge, by expanding Eq. (12) to lowest nontrivial order in , we obtain (up to a constant energy shift ) the Schrödinger Hamiltonian

 HS=∑l=1,2(p2l2m+V1b(l))+V2b, (13)

with the mass . Equation (13) describes a center in a 2D semiconductor quantum well and has been studied in Refs. phelps (); louie (); larsen92a (); larsen92b (); sandler (); ivanov (). One can therefore regard in Eq. (12) as a natural relativistic generalization of the impurity center problem.

In the remainder of this paper, we shall employ units with .

## Iii Variational approach

The Hamiltonian acts in the tensor space of two copies of the eight-dimensional single-particle Hilbert space. For the non-interacting system with , two-particle spinor wave functions for bound states are written as antisymmetrized products of single-particle wave functions,

 Φ(r1,r2)=A[Ψ(1)τ1,s1(r1)⊗Ψ(2)τ2,s2(r2)], (14)

where is the antisymmetrization operator and is a single-particle eigenstate with eigenenergy for the 2D relativistic hydrogen problem described by . To keep the paper self-contained, we summarize the well-known solution of in App. A. Eigenstates are labelled by the principal quantum number and by the half-integer angular momentum . The ground state of is realized for and . Hence the noninteracting () two-particle ground state has the energy , where both particles occupy the respective single-particle ground state. This two-particle state has finite total angular momentum , where the angular momentum operator is with .

To study the ground state of the interacting system, one could attempt to treat the two-body Coulomb repulsion by perturbation theory. However, one then finds that for , the resulting binding energy is always negative. In other words, first-order perturbation theory in the Coulomb interaction incorrectly predicts that there is no bound state (see below), and one has to proceed in a nonperturbative manner to investigate this issue. We here employ a variational treatment and construct a relativistic version of the Chandrasekhar wave function (2), dubbed Chandrasekhar-Dirac ansatz. We shall see that the corresponding energy functional is bounded from below and thus provides a variational estimate of the binding energy.

In this section, we focus on a valley- and spin-independent band gap, see Eq. (7). The topological band gap term in Eq. (8) then only requires a few adjustments, see Sec. IV.2. Moreover, we assume here that both quasiparticles are in the same valley but show in Sec. IV.3 that the variational result does not change for quasiparticles in opposite valley states. Since the Hamiltonian commutes with both and , where is (twice) the total spin operator, we have a spin singlet and a spin triplet state, where in the first (second) case, the spatial part of the wave function must be symmetric (antisymmetric). We shall see that for , the Chandrasekhar-Dirac ansatz predicts a bound state in the singlet but not in the triplet channel. However, for , bound states are found in both cases.

### iii.1 Chandrasekhar-Dirac ansatz

Our Chandrasekhar-Dirac ansatz is formulated as follows. We assume that the two-particle wave function has a factorized form, where is the normalized spin part (singlet or triplet) and the spatial part,

 Φ(r1,r2)=ΨI(r1)ΨO(r2)−ϵΨO(r1)ΨI(r2), (15)

where and are the normalized ground-state eigenspinors of the 2D relativistic hydrogen problem, see App. A, with replaced by variational parameters and , respectively. The parameter corresponds to the spin singlet/triplet sector, and the valley part of the wave function is understood. We mention in passing that in the triplet case, a wave function composed of two ground-state single-particle orbitals might not represent the optimal choice, see Sec. IV.1.

Since the single-particle Hamiltonian does not depend on the spin projection, we can use the same wave function for both particles. Explicitly, the spinors with have the form (see App. A)

 Ψλ(r)=Nλrγλ−1/2e−pλr(1iκλeiθ), (16)

where

 γλ = √14−Z2λα2,pλ=2ΔZλα, κλ = √1−2γλ1+2γλ=Zλα12+γλ=12−γλZλα. (17)

The normalization constant is given by

 Nλ= ⎷(2pλ)2γλ+12π(1+κ2λ)Γ(2γλ+1) (18)

with the Gamma function . Equation (16) represents the ground-state eigenspinor of the single-particle Dirac Hamiltonian

 Hλ=−iσ⋅∇+Δσz−Zλαr, (19)

with eigenvalue Note that the spinors are not orthogonal. Their overlap is given by

 S = (1+κIκO)√(1+κ2I)(1+κ2O)Γ(γI+γO+1)√Γ(2γI+1)Γ(2γO+1) (20) × (2pI)γI+1/2(2pO)γO+1/2(pI+pO)γI+γO+1.

### iii.2 Energy functional

We now evaluate the energy functional

 Eϵ(ZI,ZO)=⟨Φtot|H+|Φtot⟩⟨Φtot|Φtot⟩=⟨Φ|H+|Φ⟩⟨Φ|Φ⟩, (21)

with in Eq. (12). The index in refers to the spin singlet/triplet state, and the normalization factor is given by

 ⟨Φtot|Φtot⟩ = ⟨Φ|Φ⟩=2⟨ΨI|ΨI⟩⟨ΨO|ΨO⟩−2ϵ⟨ΨI|ΨO⟩2 (22) = 2(1−ϵS2).

Next, the matrix element of the single-particle Hamiltonian has the form

 ∑l=1,2⟨Φtot|HD(l)|Φtot⟩= (23) =2(∑λ⟨Ψλ|HD|Ψλ⟩−2ϵ⟨ΨI|HD|ΨO⟩S),

where we have used the normalization of the one-particle spinors. By writing the single-particle Hamiltonian as

 HD=Hλ+(Zλ−Z)αr, (24)

one can directly evaluate the matrix elements. We obtain

 ⟨Ψλ|HD|Ψλ⟩ = 2Δγλ+(Zλ−Z)Vλ, (25) ⟨Ψ¯λ|HD|Ψλ⟩ = 2ΔγλS+(Zλ−Z)U = 2Δγ¯λS+(Z¯λ−Z)U,

where for and vice versa, and

 Vλ = ⟨Ψλ|(α/r1)|Ψλ⟩=αpλγλ, (26) U = ⟨ΨI|(α/r1)|ΨO⟩=αpI+pOγI+γOS.

It is reassuring to note that the minimum with respect to for the energy

 ⟨Ψλ|HD|Ψλ⟩=2Δ√1/4−Z2λα2+2ΔZλ(Zλ−Z)α2√1/4−Z2λα2 (27)

occurs exactly at provided that . Since we have used the exact structure of the Dirac-Coulomb wave function, the result reproduces the exact ground-state energy, .

We now proceed with the two-body matrix element, see also App. B,

 V2b = ⟨Φtot|(α/r12)|Φtot⟩=2(Vdir2b−ϵVexc2b), (28) Vdir2b = ∫dr1dr2|ΨI(r1)|2αr12|ΨO(r2)|2, Vexc2b = ∫dr1dr2[Ψ†I(r1)ΨO(r1)]αr12[Ψ†O(r2)ΨI(r2)].

The standard procedure to evaluate the multiple integrals in and is to use a 2D partial wave expansion of . However, the resulting series representation converges only very slowly. Following Ref. sandler (), we found it more convenient to use an integral representation in terms of elliptic functions,

 Vdir2b = 2απp2γI+1Ip2γO+1O(γI+γO+12)B(2γI+1,2γO+1)∫10dsK(s)[s2γO(pI+spO)2(γI+γO)+1+s2γI(spI+pO)2(γI+γO)+1], (29) Vexc2b = 8απ(pI+pO)S2(γI+γO)B(γI+γO,γI+γO)∫10dsK(s)sγI+γO(1+s)2(γI+γO)+1, (30)

where is Euler’s beta function and denotes the complete elliptic integral of first kind olver ().

Importantly, we here calculated the matrix elements of the full interaction operator rather than those of the projected operator, , which are more difficult to obtain and would require a detailed numerical analysis. Both matrix elements coincide if the trial wave function has vanishing projection onto the negative energy eigenfunctions of . We show in Sec. III.3 that this is in general not the case, and hence using the unprojected Coulomb interaction is strictly speaking not justified. However, the energy functional turns out to be bounded from below and does predict a two-particle bound state. More importantly, we have verified that for the cumulative weight of negative energy states in our trial wave function is very small (), see Sec. III.3. Indeed, negative energy states will only be important if typical interaction matrix elements can overcome the band gap . For small , one therefore expects at most small quantitative corrections in the bound-state energy because of this approximation. For a treatment of stronger interactions with , however, one must resort to the matrix elements of the projected two-body operator. We leave this task for future work.

We now collect all terms and obtain

 Eϵ(ZI,ZO) = ∑λ(2Δγλ+(Zλ−Z)(Vλ−ϵSU)1−ϵS2) (31) + Vdir2b−ϵVexc2b1−ϵS2.

The energy functional has the following interesting features. First of all, is symmetric under an exchange of its arguments. Second, as illustrated in Fig. 1 for the spin-singlet case, this energy is bounded from below as long as . Third, for small , we have checked that reduces to the corresponding nonrelativistic energy functional for the problem in 2D semiconductors sandler (). However, in contrast to the nonrelativistic case, is not homogeneous in , and hence the bound-state energy explicitly depends on . As in the nonrelativistic case, this energy minimum is realized for unequal values of and .

With , the binding energy of the two-body bound state is defined for the optimized choice of as

 Eϵ,b(ZI,ZO) = Δ(1+2γ)−Eϵ(ZI,ZO) (32) ≡ α2Δ2¯Eϵ,b(ZI,ZO),

where denotes the energy of a state in which one of the particles is in the ground state of and the other in the lowest positive energy state of the continuum spectrum of , just above the gap. Equation (32) defines the rescaled dimensionless binding energy (in units of ). In the singlet case, approaches the nonrelativistic value sandler () for . For finite , deviations of from indicate the importance of relativistic effects.

Figure 1 shows that the energy functional for the singlet state with has a minimum for all studied values of . Moreover, the energy minimum is located below the threshold, i.e., the binding energy is positive, and we have a two-body bound state in the spin singlet sector. In contrast to that, our variational approach predicts that the energy functional has a minimum also for the spin triplet sector, but the minimum is now above the threshold and thus does not describe a bound state. We also notice that the energy functional for simply yields the ground-state energy of a two-particle state with the Coulomb repulsion treated within first-order perturbation theory. As anticipated above, we find , and therefore perturbation theory is not able to correctly describe the bound state for .

Table 1 lists the binding energy in the spin singlet sector for several values of . Note that both and are (well) below the critical value for all cases considered in Table 1. The fact that for the minimum happens to be at and (or vice versa, due to the symmetry of ) is rationalized by noting that only in this case the two factors and in Eq. (31) will have opposite sign. Physically, one quasiparticle then partially screens the impurity charge seen by the other quasiparticle.

We observe from Table 1 that relativistic effects tend to increase the binding energy. Approximately, we find the scaling illustrated in Fig. 2. It is interesting to note that the variational parameter increases with while decreases, see Fig. 3. Since the atomic Bohr radius is , we conclude that in the relativistic case the outer (inner) electron will be closer to (further away from) the nucleus than in the nonrelativistic case.

### iii.3 Validity of the variational approach

Before we conclude this section, we comment on the validity of this variational calculation. First, we observe that the single-particle orbital is the normalized ground state of a modified Dirac Hamiltonian which is related to by

 HD=Hλ+(Zλ−Z)αr. (33)

Assuming that is in the continuous spectrum of , we have

 (E−Eλ)⟨Ψλ|ΨE,j⟩=⟨Ψλ|[(Zλ−Z)α/r]|ΨE,j⟩. (34)

As a consequence, for and , the state is orthogonal to . This is to be expected since in that case, they are eigenstates of the same Hermitian operator with different eigenvalues. However, for , both states will generically have a finite overlap. Since we take as the ground-state orbital, it is clear that the overlap with states in the negative-energy continuum will be suppressed by a factor . Actually, the overlap can be evaluated analytically, see App. A. We can thereby compute the total weight of our variational wave function on the negative energy states,

 Wλ(α)=∫E<−ΔdE2π√E2−Δ2∣∣⟨Ψλ|ΨE,j⟩∣∣2. (35)

The result is shown in Fig. 4 for several values of and the corresponding from Table 1. We find that the total weight for is at most of order . It then stands to reason that neglecting the projection operator in the evaluation of Coulomb matrix elements does not significantly affect the variational estimate of the binding energy for .

## Iv Generalizations

So far we have restricted ourselves to the case of two quasiparticles in the same valley, in the presence of a band gap, and for an impurity of charge . In this section, we briefly address various extensions, namely (i) the case of an impurity with charge , (ii) a topological band gap, and (iii) quasiparticles in different valleys.

### iv.1 Impurity charge Z>1

From a theoretical perspective, the case is less interesting than because a bound state is then found already in perturbation theory. In fact, when treating the two-body Coulomb interaction perturbatively, the energy of the lowest singlet state, taken in the simple factorized form , coincides with the value of the energy functional for and ,

 Epert=E0(Z,Z)=4Δγ+Vdir2b. (36)

Straightforward evaluation of Eq. (36) then shows that the perturbative estimate for the binding energy , cf. Eq. (32), is negative for but positive for . Hence a bound state is predicted already by perturbation theory for , in contrast to the case .

Moreover, for , our variational method turns out to be restricted to rather small values of . Table 2 summarizes the rescaled binding energies, , and the optimal charges, , for and several values of . We here consider only the regime to ensure that the optimal charges remain well below the singular value .

For , we also find a bound state in the spin-triplet channel. However, the variational wave function used here for the triplet state is probably not the most appropriate. In the case of the helium atom, textbooks bransden () show that the simplest wave function for the lowest triplet state combines the single-particle ground state and the first excited state. By analogy, for our 2D case, a better choice for the triplet case might be to take instead of Eq. (15) the ansatz

 Φj=0=Ψ0,12,I(r1)Ψ1,−12,O(r2)−Ψ1,−12,O(r1)Ψ0,12,I(r2), (37)

which is an eigenstate of the total angular momentum operator with eigenvalue . Another option is

 Φj=1=Ψ0,12,I(r1)Ψ1,12,O(r2)−Ψ1,12,O(r1)Ψ0,12,I(r2). (38)

In the absence of Coulomb interactions, this state has the same energy as but Coulomb interactions will mix them. Therefore a variational approach should take into account both of them. However, this analysis goes beyond the scope of this paper.

### iv.2 Topological band gap

Let us next consider the case of a topological gap, , see Eq. (8), where we set as we shall still assume that both quasiparticles are in the same valley. (For related studies of the case without Coulomb impurity, see Refs. rodin (); zhong (); trushin ().) Since the total Hamiltonian now does not commute with anymore, we must distinguish whether the two quasiparticles have same or opposite spin projections, or . If both quasiparticles have the same spin projection (e.g., ), we are back to the case discussed in Sec. III but with (spin triplet). Within our variational approach for , there is no stable bound state.

Turning now to , we cannot separately consider spin singlet and triplet states. A natural way to construct the Chandrasekhar-Dirac variational wave function is as follows. We consider the 2D subspace spanned by the wave functions

 Φ1 = ΨI,↑(r1)|↑⟩1⊗ΨO,↓(r2)|↓⟩2 (39) −ΨO,↓(r1)|↓⟩1⊗ΨI,↑(r2)|↑⟩2, Φ2 = ΨI,↓(r1)|↓⟩1⊗ΨO,↑(r2)|↑⟩2 −ΨO,↑(r1)|↑⟩1⊗ΨI,↓(r2)|↓⟩2,

where spin and orbital degrees of freedom are entangled. In Eq. (39), is the eigenstate of with eigenvalue for particle , and refers to the normalized ground state of the single-particle Hamiltonian

 Hλ,s=−iσ⋅∇+sΔσz−Zλαr, (40)

with . Again, are variational parameters, and the valley part is kept implicit. Since the valley part is assumed symmetric, we have to choose antisymmetric combinations in Eq. (39).

With the ground state of given by

 Ψλ,↑∼rγλ−1/2e−pλr(1iκλeiθ), (41)

we directly obtain

 Ψλ,↓=−iσy(Ψλ,↑)∗∼rγλ−1/2e−pλr(iκλe−iθ1) (42)

as ground state of with the same energy, where . In the subspace spanned by the states , the Hamiltonian eigenvalue problem thus reduces to the problem of finding solutions of the secular equation

 det(H11−2EH12−Σ12EH21−Σ21EH22−2E)=0, (43)

where we use the notation (with )

 Hab=⟨Φa|H|Φb⟩,Σab=⟨Φa|Φb⟩. (44)

Using the results of Sec. III and noting that the single-particle matrix elements are independent of the spin projection, we find and . The roots of the secular equation are given by

 E±(ZI,ZO)=2H11−H12Σ12±|2H12−H11Σ12|4−Σ212, (45)

It turns out that this expression has the same structure as the energy functional in Eq. (31), and the corresponding results therefore apply again. We conclude that the topological band gap caused by Eq. (8) does not imply different bound-state energies as compared to the topologically trivial band gap resulting from Eq. (7).

### iv.3 Different valleys

So far we have assumed that the two quasiparticles occupy the same valley state, and we thus only have a trivial double degeneracy of the bound state. We shall now briefly discuss the case in which the two quasiparticles live in different valleys (cf. Ref. entin () for a setting without Coulomb impurity). To properly address this situation, we first recall that the Dirac-Weyl spinors represent the slowly varying parts of the electronic wave function. The complete wave function is obtained by multiplying these spinors by the appropriate Bloch wave at the or point. In the continuum description, one neglects the overlap between wave functions in opposite valleys, . (Going beyond this approximation would require a study of the lattice model.) Second, we observe that in Eq. (5) does not commute with the total squared valley operator , where . This fact must be taken into account when building the appropriate Chandrasekhar-Dirac wave function. Finally, the two-body Coulomb interaction potential has the same form as for two quasiparticles with the same valley quantum number, see App. B. A straightforward calculation as in Sec. IV.2 then shows that the resulting energy functional has the same structure as before. Therefore the optimal binding energy coincides with the one for both quasiparticles in the same valley.

## V Observables

In this section, we turn to a discussion of the probability density and of the pair distribution function for the bound state. We consider the most interesting case and focus on the Dirac band gap term in Eq. (7) for the two-body spin singlet state. However, using the results in Sec. IV, it is straightforward to obtain corresponding results also for other cases of interest.

### v.1 Probability density

We start by calculating the probability density for the two-particle bound state. The density operator is and thus the probability density is given by

 ρ(r) = 1⟨Φ|Φ⟩∫dr1dr2Φ†(r1,r2)^ρ(r)Φ(r1,r2) (46) = ∫dr2|Φ(r,r2)|21+S2.

The integral can be evaluated exactly, where the result for does not depend on the polar angle . It is then convenient to consider the radial density,

 P(r) = r∫dθρ(r)=11+S2[∑λ(2pλ)2γλ+1Γ(2γλ+1)r2γλe−2pλr (47) + 2S2(pI+pO)γI+γO+1Γ(γI+γO+1)rγI+γOe−(pI+pO)r],

with normalization . The result is illustrated for and several values of in Fig. 5. With respect to the nonrelativistic case, we observe that relativistic effects tend to enhance the probability at small distance from the Coulomb impurity. Since the radial density can be probed by STM techniques, the result can be matched to the analytical result in Eq. (47). Thereby one can hope to extract, e.g., the value of the fine structure constant .

### v.2 Pair distribution function

Next we turn to the pair distribution function of the two-particle bound state,

 g(r) = 12⟨∑i≠jδ(r−(ri−rj)⟩ (48) = ∫dr1|Φ(r1−r/2,r1+r/2)|2⟨Φ|Φ⟩.

The result again turns out to be independent of the polar angle, , and is illustrated in Fig. 6. The pair distribution function can be experimentally obtained by a statistical analysis of STM images, see, e.g., Ref. Ertl (), and can provide additional information about the existence and the properties of the two-body bound state predicted here.

## Vi Concluding remarks

In this work, we have studied the two-particle bound state problem for gapped graphene in the presence of a Coulomb impurity. We have shown that a variational approach, using the projected Hamiltonian and Chandrasekhar-Dirac spinors as trial wave functions, predicts the existence of at least one bound state. We found that in contrast to the Schrödinger case, the variational energy functional is not a homogeneous function of the coupling constant . As a consequence, the optimal values of the variational parameters depend on , and the optimal binding energy has a more complicated functional dependence on . In particular, the binding energy increases with respect to the nonrelativistic case. Moreover, we have determined the relativistic corrections to the probability density and to the pair probability density. The predicted two-body bound state can thereby be accessed experimentally, e.g., by means of STM techniques.

Finally, as a possibility for future theoretical work, it would be interesting to diagonalize the projected many-body Hamiltonian for graphene in a large basis set in order to compute the ground state energy without recourse to variational wave functions. This route can also provide information about low-lying excited resonant levels, which in turn are expected to exhibit Fano lineshapes when probed in transport or by STM methods.

###### Acknowledgements.
We thank H. Siedentop for discussions. This work was supported by the network SPP 1459 (Grant No. EG 96/8-2) of the Deutsche Forschungsgemeinschaft (Bonn).

## Appendix A Relativistic 2D hydrogen atom

In order to keep the paper self-contained, we here collect known results for the single-particle Dirac-Coulomb problem in graphene, see, e.g., Refs. rmp1 (); kotov (); novikov (); gamayun (). The massive Dirac-Weyl Hamiltonian with a Coulomb impurity of charge reads ()

 HD=−iσ⋅∇+Δσz−Zαr. (49)

The Hamiltonian (49) is exactly solvable, and the bound-state orbitals can be found, e.g., in Refs. novikov (); gamayun (). Following the notation of Ref. novikov (), in polar coordinates () they are given by

 Ψn,j(r,θ)=Nn,jργ−12e−ρ2⎛⎝(φ1+cφ2)ei(j−12)θiκ(φ1−cφ2)ei(j+12)θ⎞⎠. (50)

The half-integer index denotes the eigenvalue of the total angular momentum operator . The integer index is the principal quantum number, where the energy eigenvalues are given by

 En,j=Δ√1+Z2α2(n+γ)2. (51)

We use the notation

 ρ = 2pr,p=√Δ2−E2, (52) γ = √j2−Z2α2,κ=√Δ−EΔ+E, c = γ−ZαEpj+ZαΔp=j−ZαΔpγ+ZαEp.

The functions in Eq. (50) are confluent hypergeometric functions of the first kind olver (),

 φ1(ρ) = M(γ−ZαE/p,2γ+1,ρ), (53) φ2(ρ) = M(γ+1−ZαE/p,2γ+1,ρ).

Finally, is a normalization constant such that