# Long-range exchange interaction between magnetic impurities in graphene

###### Abstract

The effective spin exchange coupling between impurities (adatoms) on graphene mediated by conduction electrons is studied as a function of the strength of the potential part of the on-site energy of the electron-adatom interaction. With increasing , the exchange coupling becomes long-range, determined largely by the impurity levels with energies close to the Dirac points. When adatoms reside on opposite sublattices, their exchange coupling, normally antiferromagnetic, becomes ferromagnetic and resonantly enhanced at a specific distance where an impurity level crosses the Dirac point.

## I Introduction

Among possible technological promises of grapheneCN () are both electronic and magnetic applications. The former include transistors and require control of graphene’s conduction, while the latter aim to build memory devices and hinge on the ability to create and manipulate local magnetic moments. Nonetheless, both these avenues actively explore the possibilities of controlling graphene properties with chemical dopants, such as hydrogen.

In electronic applications the role of dopants is to suppress otherwise strong metallic conductivity of the materialBME (); ENM (); BJN (); DQF (). Magnetic impuritiesLFM (); PGS (); YH (); Yaz (); BDS (); SAS (); EWG (); EWT (); GEW () typically interact with the conduction band and produce electronic states that carry magnetic moments and could be spreadGGM () over many lattice spacings . Of particular interest is the mutual interaction between such impurities. The potential part of this interaction is important in determining the equilibrium spatial arrangement of the dopantsCSA (); ASL (); KCA (). At the same time, the type of collective magnetic properties is sensitive to the effective exchange coupling between dopants. Both the potential part of the effective interaction energy and its spin-dependent part are mediated by the conduction -electrons of graphene.

One particularly promising dopant is hydrogen. Because it has an energy level close to the Dirac point of the conduction -band of graphene WKL () coupling of conduction electrons has a resonant character, whose scattering amplitude resembles that of a strong substitution impurityMM (). Motivated by this similarity, we are going to concentrate on the substitution model, where the interaction of the impurity with conduction electrons,

(1) |

has both the on-site potential energy and the spin part described by the exchange coupling constant ; the spin operator of the impurity couples to the local value of the conduction electron spin density .

When two dopants reside above carbon lattice sites separated by the radius-vector , their effective electron-mediated interaction, likewise, has two parts,

(2) |

The potential part has been studied both for weakLTM () and strongSAL (); LTM () impurity strength . To the contrary, the indirect spin exchange in graphene, while extensively studied perturbatively within the usual RKKY approachSar (); BFS (); BS (); GKF (); SS (); Kog (); PF (), has not been addressed for impurities with large . Two notable previous research directions should be mentioned in this regard. In Ref. KRA, the indirect exchange between resonant Anderson impurities in graphene was studied numerically with the emphasis on the short-distance behavior. No analytic dependence has been reported, however, in the strong coupling limit. In Refs. BB, ; CZW, the indirect spin exchange in the strong limit has been addressed in a somewhat reminiscent situation of a topological insulator, but two related features distinguish graphene from that case.

First, in the low-energy Hamiltonian of a topological insulator, , the spin operator is essentially the same as (or at least directly related to) the corresponding operator in Eq. (1). In graphene, in contrast, the Pauli matrices in the Hamiltonian relate to the pseudospin operator, acting in the sublattice space. Second, in graphene the interference of the electron states from the two Dirac points results in different signs of the inter-dopant interaction when the two dopants resite on the opposite sublattices (AB-case) as compared with the same sublattice (AA) arrangement. This complication did not arise in Refs. BB, ; CZW, . In the present paper we are going to investigate the indirect spin exchange coupling between two dopants residing directly above carbon atoms of intrinsic graphene. We study how depends on 1) the sublattice arrangement of the dopants, and 2) the strength of the potential coupling , not assumed to be small. At the same time, it will be sufficient to limit our analysis to the lowest (second) order in the coupling constant , which is presumably never too large.

Let us first recapitulate the existing results for intrinsic graphene with the Fermi level at the Dirac points. To the lowest second order in and , the results are rather identical (barring the replacement ) for both and . When the two dopants reside on the same sublattice, both quantities are and are negative, meaning that the dopants attract and favor ferromagnetic arrangement of their spins. On the other hand, the sign of both interactions is reversed in the AB-case: and , indicating repulsion and anti-ferromagnetic coupling. The dependence on the distance remains the same but the overall magnitude is three times stronger than in the AA-case.

When the strength of the impurity becomes large, , the potential part changes sign, compared with the weak coupling limit. This happens for both AA and AB arrangements. This sign reversal is the feature of the linear Dirac spectrum. Additionally, the dependance of the inter-dopant interaction becomes long-rangeSAL (); LTM () (up to additional logarithmic factors). It turns out that the interaction energy in the AB-case is no longer stronger by a mere numerical factor, compared with the AA configuration (as is the case for weak ), but exceeds it logarithmically, by a large factor .

Our goal is to perform a corresponding analysis of . Below we carry it out for both sublattice configurations and for any strength of the potential assuming only that . In Section II the energy spectrum in the presence of two impurities is discussed. In Section III the general expressions for the effective exchange constant are derived in terms of the integrals over the energy of conduction electrons. In Section IV those general expressions are evaluated in the strong coupling limit.

## Ii Energy levels of a two-impurity system

To understand how and why the effective interaction of the dopants is sensitive to their relative positions on the two sublattices, it is necessary to discuss the differences in the electronic spectra induced by the presence of two dopants with a strong on-site potential coupling , elucidated in Ref. LTM, . When only one potential impurity is present, it may induce a low-energy state, provided that is strong enough. The energy of the state is determined by the equation,

(3) |

where is the area of a graphene unit cell; and the bandwidth of the conduction band . The energy if the on-site energy is large: . The energy level is quasi-localized – the overlap with the conduction band causes it to have a finite lifetime,

(4) |

Note that the impurity level is located in the lower Dirac cone, , for repulsive potentials , and vice versa.

When two impurities are present, their energy levels split. For the impurities located on the same sublattice (AA), the two energies are given by the equation,

(5) |

where the parameter is equal to the phase difference that the states belonging to different Dirac cones acquire when they travel between the two impurities. It depends on both the length of the radius-vector and the angle it makes with the zigzag direction, see Fig. 1. The oscillations in the right-hand side described by this phase, therefore, are caused by the interference of the wave functions of different Dirac species.

Importantly, both solutions of Eq. (II) remain on the same side of the Dirac point as the single-impurity level (4): no level can “escape” from its Dirac cone. This is most simply seen from the fact that no solution of Eq. (II) can ever have zero energy, .

To the contrary, when two impurities reside on different sublattices (AB configuration), the energy levels are determined from the equation,

(6) |

with the new phase being . The main new feature of the AB case is the existence of the solution at a specific distances . Such solution, from Eq. (6), has the energy (for positive ),

(7) |

and passes from the lower Dirac cone to the upper cone as the distance between the adatoms becomes shorter than a resonant distance . Right at one of the impurity levels lies exactly at the energy . Additionally, the width of this level becomes vanishingly small. As we are going to see below, the effective exchange coupling between adatoms spins is resonantly enhanced when they are separated by the distance .

## Iii General expression for the indirect exchange coupling

The Hamiltonian of two impurities (adatoms) on graphene, one positioned above a carbon atom at the origin (the sublattice it belongs to is called ) and the second above the atom some distance away (which could belong to either of the two sublattices), consists of three parts, , namely,

(8) |

Here is the kinetic energy of electrons, the hopping integral being assumed real and positive; U is the additional on-site potential energy induced by an impurity; J is the exchange coupling of an impurity spin with the spin of the conduction electrons. The hat above the electron operator indicates a spinor, the summation over the spin indices is implied; are the Pauli matrices acting in the spin (not pseudospin) space.

Because the spin exchange is never very strong, it is sufficient to determine the effective impurity-impurity spin coupling to the lowest (second) order in . Most simply this can be done by using the standard quantum-mechanical theorem that poses that the derivative of the impurity-impurity interaction energy (2) with respect to the coupling constant is equal to the expectation value of the corresponding derivative of the Hamiltonian (III):

(9) |

where we introduced Green’s function of the system, which in the interaction representation with respect to the spin part of the Hamiltonian is (spin indices shown explicitly)

(10) |

with the interaction matrix given by the standard expression,

(11) |

Expanding the -matrix to the first order in the spin Hamiltonian, we obtain,

(12) |

Here, is Green’s function of the system in the absence of the exchange coupling. Note that this function is exact with respect to the potential coupling which is not presumed to be weak. The problem of finding in the presence of two impurities was solved in Ref. LTM, . The following identity expresses it in terms of free electron () Green’s function :

(13) |

The -matrix describes the renormalization of the impurity strength from multiple scattering events,

(14) |

whereas denotes Green’s function at coinciding points,

(15) |

The area of a graphene unit cell is denoted with while the Dirac velocity is . The transposed function is given by the same expression as Eq. (13), where one simply replaces .

Integrating Eq. (12) with respect to J, and substituting the expressions for Green’s functions, we obtain the effective exchange coupling constant in the integral form,

(16) |

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

(17) |

In writing Eq. (16) we utilized that time-ordered Green’s functions do not have singularities in the first and third quadrants of the complex plane and rotated the integration path from the real axis counterclockwise to the angle so that it coincides with the imaginary axis, .

Conveniently, this removes the imaginary part of the logarithm in Eq. (15), which now becomes,

(18) |

To determine the coupling constant given by Eq. (16) it only remains to ascertain the product of two Green’s functions, , whose meaning is rather transparent: it determines the likelihood of a conduction electron to propagate from one impurity to another and then back to the first one.

The approximation (18) is valid provided that the energies are much smaller than the bandwidth, . If, in addition, the distance between the impurities is much larger than the lattice spacing, , a small-momentum expansion around the two Dirac points in the Brillouin zone is applicable to the calculation of Green’s function. The corresponding calculations have been performed previously in Refs.SAL, ; LTM, . The result is very sensitive to the positions of the two points. If both belong to the same sublattice, one obtains

(19) |

where is the Macdonald function of the zeroth order.

When two locations belong to the opposite sublattices, a different expression is found,

(20) |

where is the Macdonald function of the first order. A different sign, compared with Eq. (19), is the result of quantum interference, which is responsible for the opposite signs of interaction in AA and AB configurations.

Our general result (16) can be illustrated first by calculating the exchange coupling for the weak limit, , where the denominator in Eq. (16) can be ignored. As a result, one finds ferromagnetic coupling for the AA configuration of the two impurities,

(21) |

For the AB configuration the coupling is anti-ferromagnetic (and stronger by a numerical factor),

(22) |

The dependence of the two exchange coupling constants, as well as the signs, are consistent with the results in the existing literature, though the additional cosine and sine factors have been overlooked so far. We note that simply replacing with also yields the correct expressionsLTM () for the potential part of the impurity interaction, . This is not surprising, since the perturbation calculations for the two parts of the interaction are identical.

## Iv Strong limit

The perturbative couplings (21)-(22) are rather short-range falling off as the third power of the distance. However, it turns out that as the potential part of the electron-impurity interaction increases, so does the range of the exchange coupling. As evident from the form of the -matrix, the perturbation theory fails if , where the integral in Eq. (16) has to be calculated differently.

### iv.1 AA configuration

Let us first consider the situation of both adatoms residing on the same sublattice. Since low energies are important at large distances, , we can approximate the Macdonald function as and write with the help of Eqs. (19) and (16),

(23) |

where we introduced the dimensionless distance . In the weak coupling limit this distance is large, . In contrast, in the strong coupling limit this parameter is small. The integral in Eq. (IV.1) converges on , thus justifying the small-argument approximation used for the Macdonald function .

To calculate the integrals in Eq. (IV.1) in the logarithmic approximation, we notice that the logarithms in the integrand are both large and slow functions of their arguments. It is therefore tempting to use the standard approach to such integrals and approximate the logarithms with their fixed values taken at the characteristic arguments of the integrand, . It is easy to see, however, that this would lead to the vanishing of the integral, as all the singularities would be located in the same (lower) half-plane of the complex . It is thus necessary to calculate carefully the subleading contribution that stems from the variation of the logarithms with . This calculation is presented in the Appendix. Its result is

(24) |

where we used the shorthands, , , and . The obtained result (24) simplifies in the limit of , corresponding to strong and large distances , in which case the coupling constant becomes,

(25) |

The limit might be difficult to realize (in which case the more general formula (24) ought to be used), but it illustrates two features of . First, the coupling is antiferromagnetic: exchange interaction between adatoms reverses sign compared with the weak- limit, Eq. (21), quite reminiscent of the potential part of the adatom interactionSAL (); LTM (). Second, the coupling constant decays very weakly, logarithmically only. Needless to say, this behavior extends only to the distances , where the strong coupling limit crosses over to the weak coupling, Eq. (21), and where both expressions become of the same order of magnitude (albeit of different sign).

Fig. 2 illustrates the dependence of on the distance for different values of the potential coupling strength . For weak the coupling is ferromagnetic. With increasing antiferromagnetic coupling emerges for small distances whereas at large distances ferromagnetic coupling reemerges. With the further increase of the antiferromagnetic range extends to progressively larger distances.

### iv.2 AB configuration

When the impurities reside on different sublattices, we obtain from Eqs. (20) and (16), upon utilizing the approximation for the Macdonald function, ,

(26) |

For small , the poles of the integrand now reside on the opposite sides of the real axis. The use of the low- approximation of the Macdonald function is justified by the fact that the integral converges at . In the leading logarithmic approximation it is sufficient to take the logarithm at its typical value within the interval of convergence. One thus obtains,

(27) |

Interestingly, the coupling increases with the distance between adatoms. The sign of the coupling is antiferromagnetic, similar to Eq. (22) valid at large distances where the interaction of adatoms is perturbative.

The crossover from strong coupling limit, Eq. (27), to the weak coupling limit, Eq. (22) occurs at . In fact, a resonance takes place at , which corresponds to a localized impurity level crossing over from one Dirac cone to another, see Eq. (7). Indeed, the energies of the impurity levels are determined by the zeros of the denominator of the integrand in Eq. (16).

For small values the most important contribution into the integral in (IV.2) comes from small arguments . Keeping only the lowest order terms in in the denominator of the integrand, we write for the integral in Eq. (IV.2), while denoting , in the leading logarithmic approximation,

(28) |

In usual notations Eq. (IV.2) now reads,

(29) |

The resonant coupling (29) is ferromagnetic, in contrast to the above considered limits of long distances and short distances, both of which favor antiferromagnetic ordering of impurity spins. The origin of the resonant coupling is the existence of the zero-energy states of the two impurities at distance and the ensuing increase in the scattering of conduction electrons between the impurities. Fig. 3 illustrates the exchange coupling in the AB-configuration: the weak antiferromagnetic coupling and the stronger resonant ferromagnetic interaction.

## V Summary

The impurities (adatoms) in graphene interact via the exchange of virtual electron-hole excitations. Such interaction has the potential part as well as the effective spin exchange term. The resulting coupling strength is extremely sensitive to the strength of the on-site potential energy that the conduction electrons experience when they hop on the carbon atom located above the adatom. For weak the interaction is mediated by the band states and can be treated perturbatively. It is antiferromagnetic, , and generally stronger by a numerical factor when two adatoms reside on opposite sublattices, compared with the ferromagnetic coupling, for adatoms on the same sublattice. In both instances, when is weak, the impurity states have large energies and thus play no role when the distance between adatoms significantly exceeds the lattice spacing.

To the contrary, with increasing the impurity levels move closer to the Dirac points; as a result, the adatom-adatom interaction is mediated mostly via those levels rather than through the band states of graphene. Even though in the limit of very large the spin-spin coupling vanishes, its dependence on the distance for finite becomes highly non-trivial. In the AA-configuration of adatoms, the coupling becomes very long-range decreasing only logarithmically with the distance while being antiferromagnetic in sign – opposite to the weak coupling limit. In the AB-configuration the presence of the impurity levels results in two surprising features: the appearance of the interval of distances where increases with the distance for before undergoing a sign reversal and a resonant enhancement with the maximum at , the distance where one of the impurity levels crosses the Dirac point.

###### Acknowledgements.

We thank Dima Pesin and 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 Logarithmic Integrals

To calculate the integral in Eq.(IV.1), let us rescale the integration variable, , and introduce the shorthands , , :

(30) |

To the logarithmic approximation, we utilize the fact that the integral converges at , where . Expanding the integrand up to the first order in , we arrive at,

(31) |

The last integral has the form , where is a rational function with all its singularities located in the upper half-plane of complex : this follows from , and the fact that . Defining now a new function according to , one can use the integration by parts to obtain,

(32) |

In performing this transformation we have used that since the function does not have any singularities in the lower half-plane of . Additionally, to express the principal value integral in Eq. (A) via , we observe that

(33) |

as the integral in the left-hand side of Eq. (33) is zero for the already familiar reason: all its poles reside in the upper half-plane. From Eq. (A) we obtain that the exchange coupling constant (A) is expressed in term of the following integral of a rational function,

(34) |

The integral is straightforward and ultimately reproduces Eq. (24) of the main text.

## References

- (1) A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, Rev. Mod. Phys. 81, 109 (2009)
- (2) A. Bostwick, J. L. McChesney, K. V. Emtsev, T. Seyller, K. Horn, S. D. Kevan, and E. Rotenberg, Phys. Rev. Lett. 103, 056404 (2009).
- (3) D. C. Elias, R. R. Nair, T. M. G. Mohiuddin, S. V. Morozov, P. Blake, M. P. Halsall, A. C. Ferrari, D. W. Boukhvalov, M. I. Katsnelson, A. K. Geim, K. S. Novoselov, Science 323, 610 (2009).
- (4) R. Balog, B. Jørgensen, L. Nilsson, M. Andersen, E. Rienks, M. Bianchi, M. Fanetti, E. Løgsgaard, A. Baraldi, S. Lizzit, Z. Sljivancanin, F. Besenbacher, B. Hammer, T. G. Pedersen, P. Hofmann, and L. Hornek, Nature (London) 9, 315 (2010).
- (5) J. Ding, Zh. Qiao, W. Feng, Y. Yao, and Q. Niu, Phys. Rev. B 84, 195444 (2011).
- (6) P. O. Lehtinen, A. S. Foster, Y. C. Ma, A. V. Krasheninnikov, R. M. Nieminen, Phys. Rev. Lett. 93, 187202 (2004).
- (7) V. M. Pereira, F. Guinea, J. dos Santos, N. M. R. Peres, A. H. C. Neto, Phys. Rev. Lett. 96, 036801 (2006).
- (8) O. V. Yazyev, L. Helm, Phys. Rev. B 75, 125408 (2007).
- (9) O. V. Yazyev, Rep. Prog. Phys. 73, 056501 (2010).
- (10) V. W. Brar, R. Decker, H.-M. Solowan, Y. Wang, L. Maserati, K. T. Chan, and et. al., Nat. Phys. 7, 43 (2011).
- (11) E. J. G. Santos, A. Ayuela, D. Sanchez-Portal, New Jour. Phys. 14, 043022 (2012).
- (12) T. Eelbo, M. Waśniowska, M. Gyamfi, S. Forti, U. Starke, and R. Wiesendanger, Phys. Rev. B 87, 205443 (2013).
- (13) T. Eelbo, M. Waśniowska, P. Thakur, M. Gyamfi, B. Sachs, T. O. Wehling, S. Forti, U. Starke, C. Tieg, A. I. Lichtenstein, et al., Phys. Rev. Lett. 110, 136804 (2013).
- (14) M. Gyamfi, T. Eelbo, M. Waśniowska, and R. Wiesendanger, Phys. Rev. B 84, 113403 (2011).
- (15) H. González-Herrero, J. M. Gómez-Rodriguez, P. Mallet, M. Moaied, J. J. Palacios, C. Salgado, M. M. Ugeda, J.-Y. Veuillen, F. Ynduráin, and I. Brihuega, Science 352 437 (2016).
- (16) V. V. Cheianov, O. Syljuasen, B. L. Altshuler, and V.I. Falko, Europhys. Lett. 89, 56003 (2010).
- (17) D. A. Abanin, A. V. Shytov, and L. S. Levitov, Phys. Rev. Lett. 105, 086802 (2010).
- (18) S. Kopylov, V. Cheianov, B. L. Altshuler, and V.I. Fal’ko, Phys. Rev. B 83, 201401(R) (2011).
- (19) T. O. Wehling, M. I. Katsnelson, A. I. Lichtenstein, Chem. Phys. Lett 476, 125 (2009).
- (20) V. V. Mkhitaryan and E. G. Mishchenko, Phys. Rev. B 86, 115442 (2012).
- (21) S. LeBohec, J. Talbot, and E. G. Mishchenko, Phys. Rev. B 89, 045433 (2014).
- (22) A. V. Shytov, D. A. Abanin, and L. S. Levitov, Phys. Rev. Lett. 103, 016806 (2009).
- (23) S. Saremi, Phys. Rev. B 76, 184430 (2007).
- (24) L. Brey, H. A. Fertig, and S. Das Sarma, Phys. Rev. Lett. 99, 116802 (2007).
- (25) A. M. Black-Schaffer, Phys. Rev. B 81, 205416 (2010).
- (26) I. C. Gerber, A. V. Krasheninnikov, A. S. Foster, and R. M. Nieminen, New J. Phys. 12, 113021 (2010).
- (27) M. Sherafati and S. Satpathy, Phys. Rev. B 83, 165425 (2011).
- (28) E. Kogan, Phys. Rev. B 84 115119 (2011).
- (29) S. R. Power and M. S. Ferreira, Crystals 3, 49 (2013).
- (30) I. V. Krainov, I. V. Rozhansky, N. S. Averkiev, and E. Lähderanta, Phys. Rev. B 92, 155432 (2015).
- (31) R. R. Biswas and A. V. Balatsky, Phys Rev B, 81, 233405 (2010).
- (32) H.-R. Chang, J. Zhou, S.-X. Wang, W.-Y. Shan, and Di Xiao, Phys. Rev. B 92, 241103(R).