Potential and spin-exchange interaction between Anderson impurities in graphene

# Potential and spin-exchange interaction between Anderson impurities in graphene

M. Agarwal Department of Physics and Astronomy, University of Utah, Salt Lake City, UT 84112, USA    E. G. Mishchenko Department of Physics and Astronomy, University of Utah, Salt Lake City, UT 84112, USA
###### Abstract

The effective interaction between resonant magnetic Anderson impurities in graphene, mediated by conduction electrons, is studied as a function of the strength of the onsite energy level of the impurities and the amplitude of coupling to conduction electrons. The sign and character of the interaction depend on whether the impurities reside on the same or opposite sublattices. For the same (opposite) sublattice, the potential interaction is attractive (repulsive) in the weak coupling limit with dependence on the distance; the interaction reverses sign and becomes repulsive (attractive) in the strong coupling limit and displays behavior. The spin-exchange coupling is ferromagnetic (antiferromagnetic) at both large and small distances, but reverses sign and becomes anti-ferromagnetic (ferromagnetic) for intermediate distances. For opposite sublattices, the effective spin exchange coupling is resonantly enhanced at distances where the energy levels cross the Dirac points.

## I Introduction

Doping novel two-dimensional materials with magnetic atoms is one of the active areas of research whose ultimate objective is to design systems with the desired magnetic properties. To better exploit an emerging magnetism in such doped materials, it is important to understand how magnetic impurities interact with each other.

Impurities in conventional three-dimensional metals induce famous charge (Friedel) and spin density (RKKY) oscillations of the conduction electron density, , in the long distance limit. In conventional two-dimensional electronic systemsFK (); BM () the amplitude of these oscillations decays inversely proportional to the square of the distance. One exception is graphene, a two dimensional material known for its remarkable electronic properties and a potential for applicationsCN (). Density oscillations in graphene in both intrinsic (undoped) and extrinsic (doped) limit decay asMM () , much like in three-dimensional systems. The RKKY interaction between magnetic impurities in graphene has also been studied extensivelySar (); BFS (); BS (); GKF (); SS (); Kog (); PF (); AFS (). The sign of the RKKY interaction for a bipartite lattice of intrinsic graphene at half-filling is dictated by the particle-hole symmetry and is anti-ferromagnetic (ferromagnetic) when the impurities reside on different (same) sublattices. This is found at all length scalesSar (). For example, RKKY exchange coupling between spins of impurities located on the same sublattice has the following oscillatory behaviorSS () , where and are the positions of two Dirac valleys in the reciprocal lattice. The coupling between spins on different sublattices, , has a similar oscillatory pattern, but the negative sign and the amplitude that is three times larger than in the AA case.

The above referenced studies considered interactions of impurity atoms with the host material perturbatively. On the other hand, some adatoms (such as hydrogen) are better described by the Anderson model of a localized orbital hybridized with a conduction band of a host material. Such a model allows for a strong coupling of the localized orbital to conduction electrons. In the present paper we consider two Anderson impurities with a low energy orbital hybridized with the -band of graphene with some amplitude . We further assume that the orbital is below the Fermi energy (Dirac point) of undoped graphene, , and that the Coulomb onsite energy is large enough, , so that there is always an uncompensated spin associated with the impurity. Such assumptions work well for hydrogen adatoms, which have energy level close to the Dirac point of grapheneWKL (). It is known that chemisorption of hydrogen atoms on graphene can indeed induce magnetic moments GGM ().

Magnetic applications of graphene would benefit from the ability to control magnetic moments. This, in turn, requires the knowledge of the magnitude and sign of the effective exchange coupling between dopants. Of particular interest is the behavior of resonant Anderson impurities, where the orbital resides close to the Dirac pointsMM (); WKL (). This results in the enhanced scattering of conduction electrons off the impurityMM ().

It is instructive to begin our analysis of the Anderson-type impurities with a discussion of potential impurities. Recent studies of impurity-impurity interaction in the case of substitution impurities with an onsite potential have obtained an analytical expression exact in SAL (); LTM (); AM (). In particular, in the large limit, interaction between impurities on the same sublattice is long range, (up to logarithmic factors), and is repulsive, in contrast to the weak limit where it decays as and is attractive. The interaction between impurities residing on opposite sublattices similarly reverses sign and changes behavior when the strength varies. Effectively, the Anderson impurity maps on the potential impurity model if one replaces the onsite potential strength with the energy-dependent coupling parameter , . The weak potential impurity limit, analogous to the Anderson model with small , maps onto the case of a large onsite energy and/or small amplitude such that is an energy independent constant for most energies except . As a result, the interaction of both types of impurities depends on the coordinate in a similar fashion, . With decreasing the distance between the impurities, the strong coupling limit is achieved when the on site energy becomes of the order of . In this strong coupling limit, the effective interaction energy is given by this very ratio , since once drops out of the picture, there is only one remaining low-energy scale in the system. The sign of the interaction is repulsive (attractive) for impurities belonging to the same (different) sublattices. Below, in Section III we confirm that Anderson impurities in the resonant coupling regime demonstrate a similar -dependence.

We then investigate the effective spin-spin exchange coupling between two Anderson impurities in graphene and compare it with our recent results for substitutional potential impuritiesAM (). We limit our analysis to the lowest (second) order in the coupling between the localized spins and conduction electrons but explore a broad range of the parameters and . The case of weak Anderson impurities yields, , similar to effective spin coupling in the potential impurity case. We then explore how behaves in the strong coupling limit. In our recent workAM (), we have shown that the effective spin exchange coupling between substitutional magnetic impurities can become resonantly enhanced at a specific distance where an impurity level crosses the Dirac point. We find similar enhancement for Anderson impurities for sufficiently large couplings .

This paper is organized as follows: In Section II, we discuss the energy levels of two Anderson impurities. In Section III and IV, we derive general expressions for the potential interaction and the spin exchange coupling between the impurities respectively and consider the limiting case of .

## Ii Energy levels of a two-impurity system

We consider the tight binding model of electrons in graphene interacting with two Anderson impurities located at and . In order to calculate the interaction between the impurities, we first determine the energy spectrum of the system. The Hamiltonian of the system consists of the kinetic energy of electrons, the on site impurity energy level , and a coupling term describing hybridization of conduction band with the impurity states with amplitude ,

 Ha= t∑r∑i=1,2,3^ψ†(r)^ψ(r+ai)+ϵo∑j=1,2^d†(rj)^d(rj) +γ∑j=1,2(^d†(rj)^ψ(rj)+^ψ†(rj)^d(rj)). (1)

Here is the hopping integral, is electron operator of conduction electrons; when belongs to sublattice A, when it belongs to sublattice B, and is the operator of the localized impurity states. The index enumerates the impurities. Using the Fourier representation for the electron operators, , we obtain from the equation of motion, , the following system of coupled equations (for the AB impurity configuration),

 E^a(k) = t(k)^b(k)+√2Nγd^(0), (2) E^b(k) = t∗(k)^a(k)+√2Nγe−ik⋅R^d(R), (3) E^d(0) = ϵo^d(0)+√2Nγ∑k^a(k), (4) E^d(R) = ϵo^d(R)+√2Nγ∑k^b(k)eik⋅R, (5)

where and is the total number of carbon atoms. Eliminating and gives,

 (E−ϵo)[E^a(k)−t(k)^b(k)] = 2γ2N∑k′^a(k′), (E−ϵo)[E^b(k)−t∗(k)^a(k)] = 2γ2N∑k′^b(k′)ei(k′−k)⋅R.

Solving the above two equations yields the equation for the localized impurity energy levels:

 [1−γ2∑kA(k,0)]2=γ4∑kB(k,R)∑k′B(−k′,R), (6)

where

 {A(k,R)B(k,R)}=2e−ik⋅RN(E2−|t(k)|2)(E−ϵo){Et(k)}.

The poles in the above expression should be avoided in the usual way by replacing .

Considering only low energy physics of the Hamiltonian, we expand momentum vector near the two Dirac points, . Summation over momentum vectors then gives,

 E=ϵo±α0v|sin(θAB)|/Rα0(ln|t/ϵo|+iπ2)+1,α0=γ2Aoπv2. (7)

Here , where is the angle measured from zig-zag direction as shown in Fig. 1. The dimensionless constant describes the strength of hybridization relative to the hopping integral. Importantly, one of the impurity levels in AB configuration can cross the Dirac point at a particular distance . As we will see in Sec. IV, the spin exchange coupling between impurities residing on different sublattices can become resonantly enhanced at this distance where crossing occurs.

## Iii Interaction energy: potential part

The interaction energy of conduction electrons described by the Anderson Hamiltonian (II) can be calculated using the following well-known quantum-mechanical identity,

 ∂W∂γ=⟨∂H∂γ⟩=∑j=1,2⟨^d†(rj)^ψ(rj)+^ψ†(rj)^d(rj)⟩. (8)

This identity can be written in terms of electrons Green’s function

 G(r,r′,t)=−i⟨T^ψ(r,t)^ψ†(r′,0)⟩. (9)

Because according to Eqs. (4) and (5), , we obtain from Eq. (8),

 ∂W∂γ=−2iγE−ϵo∑j=1,2G(rj,rj,t=0−). (10)

The equation for Green’s function in the energy representation is

 EGE(r,r′)−t∑iGE(r+ai,r′)−γ(E)δr,0GE(0,r′) −γ(E)δr,RGE(R,r′)=δr,r′, (11)

where we introduced the shorthand,

 γ(E)=γ2E−ϵo. (12)

The solution of Eq. (III) has been found elsewhereLTM () for the case of two impurities with the onsite potential U. Because the present case differs from that situation only by the replacement , we can use the result of Ref. LTM, ,

 GE(R,0)=GE(R,0)1+2TEGE(0,0)+T2EG2E(0,0)1−T2EG2E(0,R)GE(R,0), (13)

which expresses the two impurity Green’s function (for the electron propagation between two impurities) via the free electron Green’s function .

The interaction energy (that part of that depends on the distance between impurities) then follows from Eq. (10):

 W(R)=2i∞∫−∞dE2πln(1−T2EGE(R,0)GE(0,R)), (14)

where stands for the T-matrix,

 TE=γ(E)1−γ(E)GE(0,0). (15)

The free electron Green’s function evaluated at coinciding points, is

 GE(0,0)=−EA0πv2[ln(t|E|)+iπ2 sgnE]. (16)

Using now the fact that the time-ordered Green’s functions do not have singularities in the first and third quadrants of the complex -plane, and rotating the integration path counterclockwise by the angle so that it follows the imaginary axis, , we obtain the expression for the interaction energy,

 W(R)=−2∞∫−∞dω2πln(1−T2iωΠiω(R)), (17)

where we introduced the following shorthand for the product of two Green’s functions,

 Πiω(R)=Giω(0,R)Giω(R,0). (18)

For the AA configuration of adatomsSAL (); LTM (),

 (19)

where is the Macdonald function of the zeroth order and . For AB configuration the product is given bySAL (); LTM (),

 Πiω(R)=ω2A2oπ2v4K21(|ω|Rv)sin2θAB, (20)

where is the Macdonald function of the first order and is defined after Eq. (7).

To make subsequent calculations of the energy given by Eq. 17 more compact, let us introduce the two dimensionless parameters

 α=α01+α0ln(Ra),β=Rvϵo1+α0ln(Ra), (21)

namely, the renormalized impurity coupling strength and the parameter that characterizes the location of the impurity level relative to the energy scale of the electron travel between the impurities. With increasing the “bare” coupling , the renormalized approaches a (distance-dependent) constant. Note that in the long-range limit , to which the present theory only applies, is always less than . The parameter , as we are about to see, describes the effective strength of the impurity with large corresponding to the weak impurity limit and to the strong coupling domain.

The potential interaction energy expression for impurities residing on different sublattices in terms of and is given by,

 WAB(R)=−2vR∞∫−∞dx2πln(1−α2x2K21(|x|)sin2θAB(ix−β)2). (22)

In deriving the above expression we have used Eq. (20) for the product of two Green’s function and the expression for the matrix given by Eqs. (15) and (12). The integral in Eq. (22) can be calculated in different limits of and .

(i) When the distance between the impurities is large enough so that , we can neglect in the denominator and calculate the remaining integral by expanding the logarithms over a small ratio (as explained before, ),

 WAB(R) ≈−2vR∞∫−∞dx2πln(1−(αβ)2x2K21(|x|)sin2θAB) ≈3πα2016v3R3ϵ2osin2θAB. (23)

This is simply the weak impurity limit, where interaction decays with the distance as , just like in a case of a substitution impurity. The interaction is positive (repulsive) there. This regime is also realized when the impurity level is sufficiently far away from the Dirac point.

(ii) As the distance decreases or, alternatively, the energy level approaches the Dirac point, a situation of small is eventually realized. (This condition means that the impurity level has the energy, , i.e. much smaller than the energy corresponding to the distance ). In the most interesting case of , two scenarios can occur, depending on how compares with . In the limit of , one can still expand the logarithm in the integrand,

 WAB(R)≈vα2sin2θABπR∞∫−∞dxx2K21(|x|)(ix−β)2, (24)

even though it is no longer possible to neglect in comparison to in the denominator. Because this limit requires small , the difference between two coupling constants becomes insignificant, , whereas . The remaining integral is calculated in Appendix A to give,

 (25)

The sign of the interaction remains the same as in Eq. (III) but the dependence on changes from to a long range . As should be, the two expressions, Eq. (III) and Eq. (25), becomes of the same order when ; this happens when .

In the second limit of , where the impurity level energy is negligible in comparison with , it is appropriate to ignore in the integrand in Eq. (22),

 WAB(R) ≈−2vR∞∫0dxπln[1+α2sin2θABK21(x)] ≈−2vR∞∫0dxπln(1+α2sin2θABx2). (26)

Since the relevant values of in this integral are of the order of , neglecting in the original integral is justified. Utilizing also that is small, we used the small-argument expansion of the Macdonald function, . The remaining integral is straightforward to calculate using integration by parts and is equal to . The expression of interaction energy is thus given by,

 WAB(R)=−2αvR|sinθAB|. (27)

Note that the presence of logarithmic term in indicates the onset of multiple scattering of electrons off the impurities. This is the strong impurity limit where the interaction is attractive, in contrast to the weak impurity limit. In the limit of , we recover the expression found earlier in Refs. SAL, ; LTM, for strong potential impurities.

Fig. 2 illustrates the dependence of on the distance between impurities for different values of onsite energy : 0.1, 0.03 and 0.01 eV and coupling = 1 eV. The inset plot in the figure is to explain the behavior of interaction energy with the help of impurity strength parameter and renormalized impurity coupling . It shows the variation of with distance for the above set of onsite energies along with plotted for eV. For eV, remains greater than for all values of , hence the interaction energy is always repulsive. As we decrease to 0.03 eV and further down to 0.01 eV, we see the transition from weak coupling to strong coupling occurs leading to attractive interaction. It happens at the value of when . Fig. 3 on the other hand shows the dependence of on the distance in different regime. Even though the potential interaction in this regime is repulsive, it changes from to dependence at .

The interaction energy for the two impurities residing on same sublattices is given by,

 WAA(R)=−2vR∞∫−∞dx2πln(1+α2x2K20(|x|)cos2θAA(ix−β)2). (28)

As in the AB case, we calculate the above integral in different limits of and . In the weak impurity limit, , the integrand is simplified by neglecting in the denominator and the remaining integral can be calculated by expansion of log,

 WAA(R) ≈−πα2o16v3R3ϵ2ocos2θAA. (29)

Because the integral converges over , neglecting in comparison to in the denominator is justified. The interaction in the AA configuration in weak impurity limit is attractive, in contrast to the repulsive interaction in the AB case and is three times smaller in magnitude.

When , similarly to the AB case, two limits arise. (i) For , we are still justified to expand the logarithm (but not to neglect in the denominator),

 WAA(R)≈vα2cos2θAAπR∞∫−∞dxx2K20(|x|)(x+iβ)2. (30)

The above integral is calculated in Appendix A and the expression of interaction energy is given by,

 WAA(R)= vα2cos2θAAπR(π2−π22ϵoRv− (31)

(ii) In the remaining limit of , we can neglect in denominator of the integral in Eq. (30) and obtain,

 WAA(R)=πα2v2Rcos2θAA. (32)

We see that in the strong impurity limit the interaction is repulsive, in contrast to the weak impurity limit Eq. (III) where it is attractive.

Fig. 4 illustrates the dependence of on the distance between impurities for different values of the onsite energy : 1, 0.2, 0.01 eV and coupling eV. The inset plot in the figure is to explain the behavior of interaction energy with the help of impurity strength parameter and renormalized impurity coupling . It shows the variation of with distance for the above set of onsite energies along with plotted for eV. For eV, remains greater than one for all values of , hence the interaction energy is always attractive. As we decrease to 0.2 eV and further down to 0.01 eV, we see the transition from weak coupling to strong coupling occurs leading to repulsive interaction. It happens at the value of when . Fig. 5, on the other hand shows the dependence of on the distance in different regime. The potential interaction in this regime changes from attractive, , to repulsive, , at .

## Iv Interaction energy: spin-dependent part

To describe a spin-dependent part of the interaction between Anderson magnetic impurities in graphene, we add to our Hamiltonian the spin term,

 H=Ha+Hsp, (33)

where is given by Eq. (II) and

 Hsp=JS1⋅^ψ†α(0)^σαβ^ψβ(0)+JS2⋅^ψ†α(R)^σαβ^ψβ(R). (34)

contains two short-range exchange interactions between spins of impurities, , , and those of conduction electrons described by the Pauli matrices . The exchange coupling is assumed to be small compared with both and . As a result of the coupling to conduction electrons, there appears the effective coupling of impurity spins,

 Heff=Jeff(R)S1⋅S2. (35)

The effective exchange constant can be obtained from the already familiar method of differentiation with respect to the coupling parameter ,

 ∂Jeff∂JS1⋅S2=⟨∂H∂J⟩=∑j=1,2Sj⋅⟨^ψ†(rj)^σ^ψ(rj)⟩. (36)

The expectation values in Eq. (36) should be calculated to the lowest (first) order in the Hamiltonian (34). This yields,

 Jeff=−2iJ2∞∫−∞dE2πℏGE(R,0)GE(0,R). (37)

The last expression can be simplified further by expressing Green’s functions via free electron Green’s functions, Eq. (13). At last, rotating the integration path counterclockwise by the angle , , we obtain,

 Jeff=J2πℏ∞∫−∞dωΠiω(R)[(1−γ(iω)Giω(0,0))2−γ2(iω)Πiω(R)]2, (38)

Having obtained the general expression for the effective spin exchange coupling for two impurities in AA or AB configuration, we can now proceed with evaluating it for different limits of the coupling strength and the position of the impurity level .

### iv.1 AB configuration

We begin by calculating spin-exchange coupling in AB configuration using Eq. (20) for the product of two Green’s functions and . Green’s function at coinciding points is given by Eq. (16). To simplify calculations further let us introduce dimensionless distance as (here, , defined earlier in Eq. (7)),

 JABeff(R)=J2π3ℏα40A2ovR3sin2θAB× ∞∫−∞dxx2K21(|x|)(ix−α0ρ)4[(ρ−ix(1α0+ln[Ra|x|]))2−sin2θABx2K21(|x|)]2. (39)

Below we consider the dependence of the effective exchange constant on the distance (represented by the dimensionless variable ), as predicted by Eq. (IV.1).

At both large and small distances , we obtain the same power-law dependence, . Indeed, for very large , the integral in Eq. (IV.1) is equal to , so that the exchange coupling becomes,

 JABeff(R)=3J216πℏA2ovR3sin2θAB. (40)

The applicability of this expression follows from the observation that contribute to the integral and that must be large enough compared with , or in terms of the actual distance, .

On the other hand, at small distances one can simply set in the integral. This again reveals the dependence as the integral still converges at . Provided that , one can also neglect the term in the denominator. The exchange constant then assumes the form,

 JABeff(R)=3J216πℏA2ovR3sin2θAB(1+α0lnRa)4. (41)

The short-distance value (41) is suppressed, compared with the long-distance asymptotic (40), by the additional factor that depends weakly (logarithmically) on . To justify neglecting in both the numerator and the denominator, it is sufficient to have or, equivalently, .

At , the contribution of to the integral in Eq. (IV.1) is the only one that matters, and leads to Eq. (41). However, as increases (but still does not exceed , see below), a contribution of small might become dominant, where the term in the numerator of the integrand can be neglected in comparison with and where one can approximate . As long as the poles of the integrand are on opposite sides of the real axis (see discussion below), the remaining integral can be calculated by the residue method,

 JABeff(R)= J2π3ℏA2ovR3(ϵoRvα0)4sin2θAB× ∞∫−∞dx[(ρ−ixln[cRa|x|])2−sin2θAB]2 =π2J22π2ℏϵ4ov3RA2oγ81|sinθAB|ln(α0vaϵo). (42)

Because , the typical value and the condition is satisfied automatically. But to ensure that is smaller than , one must also have or, equivalently, . We conclude that the exchange coupling constant is the sum of contributions (41) and (IV.1) in the range of distances, . Within this range, with increasing , the contribution gradually becomes dominated by the linear term. The crossover occurs at the point where the two terms are of the same order of magnitude, at .

One additional condition should be emphasized. The residue method applies only if or equivalently ; otherwise the poles of the integrand in Eq. (IV.1) reside on the same side of the real axis. To ensure that one nonetheless has , it is necessary that the value is not too small, i.e. that .

The linear increase of the effective spin exchange coupling is caused by multiple scattering of conduction electrons off the two impurities. As the distance increases even further and approaches , a resonant enhancement of the exchange constant occurs. There, one of the energy levels of the impurities in AB configuration crosses the Dirac point, see Eq. (7). As a result, at the integrand has the singularity at . For small values the integral can be calculated by keeping only terms linear in in the denominator. This is justified by the fact that the integral converges at . In the leading logarithm approximation we obtain (see Appendix B for details),

 JABeff(R)=−vJ2ϵ2o4ℏγ4|R−R0||sinθAB|ln2(R20a(R−R0)). (43)

Fig. 6 illustrates the dependence of the spin exchange coupling on the distance between the impurities for different values of the onsite energy and eV. The interaction changes from the weak impurity limit to strong impurity limit at short distances on decreasing value of . This can be understood from the inset plot shown inside the graph, it shows the variation of and with . is always greater than for eV and thus the interaction is always anti-ferromagnetic. It changes to strong impurity limit for for eV where the impurities are ferromagnetically coupled.

Fig. 7 illustrates the dependence of the spin exchange coupling on the distance between the impurities for two values of onsite energy and eV and eV. At small distances , the interaction is anti-ferromagnetic and decreases with increasing . On increasing distance, we see a transition from weak coupling to strong coupling. Above , the interaction remains anti-ferromagnetic but increases linearly with . At , the spin exchange coupling becomes resonant ferromagnetic.

### iv.2 AA configuration

In the AA configuration of impurities the spin exchange coupling is given by,

 JAAeff(R)=−J2π3ℏα40A2ovR3cos2θAA× ∞∫−∞dxx2K20(|x|)(ix−α0ρ)4[(ρ+ix(1α0+ln[Ra|x|]))2+cos2θAAx2K20(|x|)]2. (44)

Similarly to AB case, we obtain the same dependence of spin exchange coupling for both small and large distances . For large , we can keep only the terms in both numerator and denominator and utilize the fact that the integral in Eq. (IV.2) converges at . Because , the exchange coupling becomes,

 JAAeff(R)=−J216πℏA2ovR3cos2θAA. (45)

On the other hand, for small and large value of , we can neglect in the denominator both and the term containing :

 JAAeff(R)=−J216πℏA2ovR3cos2θAA(1+α0lnRa)4. (46)

Again values determine the integral, and thus the last expression is valid as long as , which is equivalent to .

As increases, the contribution of small can become important, where we can neglect term in the numerator in comparison to . For , the integral in Eq.(IV.2) converges on where we can approximate the Macdonald function . In the leading logarithmic approximation it is typically sufficient to take the logarithms at the characteristic arguments; in this case . But the integral thus evaluated will vanish because both poles of the integrand will lie on the same side of the -axis. Thus, it is important to calculate the subleading contribution to the integral where one can no longer treat logarithms as constant. The corresponding calculation is presented in Appendix C. In the limit of small and large distances it amounts to,

 JAAeff(R)=2J2πv2ϵ3o3ℏA0γ6ln(vα0Rϵo)cos2θAAln3(vα0aϵo). (47)

We note that since , the condition requires that the value of is large, . If this is the case, the exchange coupling constant for is the sum of two contributions, Eq. (46) and Eq. (47). With increasing , the logarithmic contribution (47) exceeds the part, Eq. (46). This crossover occurs when the two terms becomes of the same order, at .

Fig. 8 illustrates dependence of the spin exchange coupling on the distance between impurities for eV and different values of the onsite energy . The impurities are ferromagnetically coupled for large distances . The change from weak to strong impurity occurs at distances, , as seen from the inset plot. The impurities are antiferromagnetically coupled in the strong impurity limit until where it switches back to ferromagnetic coupling.

## V Summary

Two singly-occupied Anderson impurities with the energy level below the Fermi energy interact resonantly with each other. The interaction is facilitated by the exchange of virtual electron-hole excitations. The sign and nature of the interaction depend on whether the impurities reside on the same sublattice or then opposite sublattices.

For opposite sublattices both the potential part of the interaction and the effective spin exchange coupling have resonant character when one of the energy levels of the two-impurity system passes through the Dirac points. The potential interaction is repulsive and decays with the third power of the distance in the weak coupling limit. The resonant potential interaction decays as the first power of the distance and is attractive. The spin-exchange part of the interaction is anti-ferromagnetic both at small and large distances. At the distances where level crosses the Dirac points, the coupling is ferromagnetic and resonantly-enhanced.

For the same sublattice, the potential part of the interaction is attractive in the weak coupling limit and repulsive in the strong coupling limit. The spin-exchange coupling is ferromagnetic at large and small distances but reverses sign and becomes anti-ferromagnetic for intermediate distances.

## Vi Acknowledgements

We thank Oleg Starykh for helpful discussions. The work was supported by the Department of Energy, Office of Basic Energy Sciences, Grant No. DE- FG02-06ER46313.

## Appendix A Calculation of integrals involving special functions

(i) The integral in Eq. (24) is of the form:

 I(β)=−∞∫−∞dxx2K21(|x|)(x+iβ)2=−2∞∫0dxx2K21(x)(x2−β2)(x2+β2)2, (48)

Integrating by parts and separating the leading contribution to the integral then gives,

 I(β)=π22+4β2∞∫0dxK1(x)(x2+β2)ddx(xK1(x)). (49)

The main contribution to the remaining integral comes from where we can expand upto second order and differentiate it to rewrite the integral as,

 I(β)= π22+4β2∞∫0dxln(0.89x)(x2+β2) =π22+2πβln(0.89β). (50)

(ii) The integral in Eq. (30) can also be calculated in a similar way. On integration by parts we get,

 I(β)=∞∫−∞dxx2K20(|x|)(x+iβ)2≈4∞∫0dxx2K0(x)(x2+β2)ddx(xK0(