Indirect Exchange and Ruderman-Kittel-Kasuya-Yosida (RKKY) Interactions in Magnetically-Doped Graphene

Magnetically-doped graphene systems are potential candidates for application in future spintronic devices. A key step is to understand the pairwise interactions between magnetic impurities embedded in graphene that are mediated by the graphene conduction electrons. A large number of studies have been undertaken to investigate the indirect exchange, or RKKY, interactions in graphene. Many of these studies report a decay rate faster than expected for a 2-dimensional material and the absence of the usual distance dependent oscillations. In this review we summarize the techniques used to calculate the interaction and present the key results obtained to date. The effects of more detailed parameterisations of the magnetic impurities and graphene host are considered, as are results obtained from ab initio calculations. Since the fast decay of the interaction presents an obstacle to spintronic applications, we focus in particular on the possibility of augmenting the interaction range by a number of methods including doping, spin precession and the application of strain.


Introduction
Graphene has been in the scientific limelight for the past few years due to a range of interesting properties that it has been shown, or predicted, to display. Superlative mechanical properties and a unique set of potentially tunable electronic and optical properties suggest a large number of possible applications for graphene and related materials [1][2][3][4]. Apart from potential roles in nanoelectronics, optoelectronics and as reinforcements in composites, graphene-based materials have also been mooted for application in spintronics [5]. This field is projected to play a very important role in future technologies as it may allow information storage, processing and communication at faster speeds, and with lower energy consumption, than is possible with conventional electronics [6][7][8][9]. The principal idea of spintronics is to exploit not only the charge of an electron, but also the spin degree of freedom associated with its intrinsic angular momentum. Despite carbon not being magnetic, graphene-based spintronics may be achievable when we consider that many of its derivative materials and nanostructures display various scenarios of magnetism [5]. Graphene has also many properties that suggest its possible use as a carrier of spin information, including weak spin-orbit and hyperfine couplings, which are the main sources of relaxation and decoherence of electron spin [10][11][12][13][14][15][16].
Many of the proposed graphene-based spintronic devices are underpinned by the spin-polarised edge state that is predicted to occur when a graphene sheet is cut to have the so-called zigzag edge geometry. Particular focus has been paid to zigzag-edged graphene nanoribbons, narrow strips of graphene with parallel zigzag edges that are predicted to display opposite spin orientations. Such a system may allow a possibility of tuning its spin-transport properties [17][18][19][20][21][22]. For example, the prospect of triggering a half-metallic state using external electric fields in zigzag-edged nanoribbons has been suggested [19]. The realisation of such a device would allow efficient electronic control of spin transport, a very useful property in spintronics and something that is difficult to achieve in other materials. Despite theoretical advances in the study of nanoribbons, and some recent experimental results hinting at signatures of this edge magnetism [23], the difficulty in patterning the edge geometries required for these effects to be observed may prevent a wider scale exploitation. Furthermore, the spin-polarised edge state for zigzag edges is predicted to be highly dependent on the edge geometry and not particularly robust under the introduction of edge disorder [24]. These factors present major obstacles in the path of utilising the intrinsic magnetic edge states of graphene in experimentally realisable devices.
Another possibility that has been proposed for graphene-based spintronics is the exploitation of defect-driven magnetic moments that arise in graphene [25][26][27]. Magnetic moments have been predicted to form around vacancies and other defects in the graphene lattice and the possibility of engineering a ferromagnetic state in graphene from such moments has been suggested. However, such a claim would seem to be restricted by the implications of the Lieb theorem [28], which states that any such magnetic moments arise from a disparity between the two sublattices of graphene. Large-scale, randomised disorder would tend to minimise such a disparity and prevent the formation of a ferromagnetic state. However, recent experimental evidence suggests the possibly of engineering such a state through partial hydrogenation [29]. The existence of such a state may then be accessible through magnetoresistance measurements [30].
A third possibility for incorporating graphene in spintronic devices, and one that we shall focus on in this review, is through the doping of graphene with magnetic impurity atoms, as shown in Figure 1. This approach allows graphene to act as host to an ensemble of transition metal atoms, or indeed more generic magnetic species, and mediate the interactions between them, potentially allowing for long ranged ordering or the transfer of spin information between impurities. In order to predict the magnetic properties of such systems, it is essential to understand both how a single impurity connects to the graphene lattice and also the nature of the interactions between multiple impurities. To address the former, a large number of studies have investigated the hybridisation and resultant electronic, magnetic and magnetotransport properties of graphene systems with different magnetic species embedded in them [31][32][33][34][35]. However, the interaction between the localised moments formed at the impurity sites is predicted to be largely independent of the magnetic impurity species chosen and instead depend strongly on the electronic structure of the host material. This is because the indirect exchange coupling (IEC) [36][37][38][39], often referred to as the Ruderman-Kittel-Kasuya-Yosida (RKKY) interaction [40][41][42][43][44], between magnetic impurities in a graphene system is mediated by the conduction electrons of the graphene host. Throughout this work we will use the terms "IEC" and "RKKY" interchangeably to refer to a general conduction-electron mediated interaction between magnetic objects, unless specifically referring to the RKKY approximation to the interaction. A long-ranged interaction of this type allows impurity moments on graphene to feel each other's presence and respond to magnetic perturbations or excitations at other impurity sites. The usual manifestation of this interaction is through the lowering of the total energy of the system when the moments adopt certain favourable alignments. In this way the interaction can dictate any magnetic ordering that arises between the moments. These types of interactions have been investigated previously in multilayer devices, leading to the discovery of the Giant Magnetoresistance effect [45,46], and more recently in materials closely related to graphene, namely carbon nanotubes [47][48][49].
To exploit the interaction between impurities, it is important to understand its behaviour as the separation between them is varied. Previous studies suggest that this type of interaction should oscillate and decay as the impurity separation is increased [36][37][38][39][40][41][42][43][44]. Furthermore, the oscillation period is determined by the Fermi surface of the host material, whereas the decay rate depends on the dimensionality of the host and magnetic objects. For impurity atoms in a one-dimensional or two-dimensional host medium, the interaction should be expected to decay as 1 D or 1 D 2 respectively, where D is the distance between the magnetic impurities. The decay rate in one-dimensional systems is borne out by investigations in metallic carbon nanotubes where the predicted long range decay is found theoretically for substitutional magnetic impurities [47]. However, unique features arising from the peculiar electronic structure of the underlying graphene lattice were also reported. The sign of the interaction, determining the preferred alignment of the moments on two impurities, was found not to oscillate as a function of distance as seen in other materials, but rather to depend only on the connections of the impurities to the two triangular sublattices composing graphene. These sublattices are represented by filled and hollow circles in Figure 1. Furthermore, for certain impurity configurations, the interaction was found to decay much faster than expected [48][49][50].
Similar features to those reported for carbon nanotubes have now been predicted in graphene. Due to the interest in carbon-based spintronics, a large number of studies in recent years have investigated indirect exchange interactions between localised moments in graphene. Sublattice-dependent signs and faster than expected decay rates have been seen when a variety of methods are employed, ranging from analytic or numerical treatments of the interaction using tight-binding or Dirac Hamiltonians to more intensive total energy calculations within a density functional theory (DFT) framework. In this review the details of the various approaches will be discussed in Section 2 and the consensus behaviour for simple impurities that has emerged from the majority of the studies analysed in Section 3. The results for impurities beyond the simplest substitutional case and for considerations often neglected by the basic RKKY approximation are discussed in Section 4. Finally, in Section 5 we consider how the interaction between magnetic impurities can be manipulated in order to extend the interaction range or otherwise alter its behaviour.

Methods
In this section we will outline several of the general methodologies employed to calculate the indirect exchange coupling in graphene. We begin by considering an energy difference calculation between the parallel and antiparallel orientations of the localised moments on two magnetic impurities located at sites A and B in the graphene lattice. These configurations are represented schematically in Figure 2.
Employing the Lloyd formula [51] to calculate the total energy difference between these two spin configurations we arrive at the Quantum Well method result proposed by Edwards et al. [38,39]. We shall briefly derive this result in a completely generic form in terms of Green functions without stating the form of the Hamiltonian chosen to describe the electronic structure of the graphene host or the magnetic impurities. From this result we identify the common approximations often made in the literature and demonstrate how the commonly used RKKY expression can be reached.
From these general results we turn our attention to the specific case of graphene. Two common Hamiltonians used to describe the electronic structure, namely the tight-binding and linearised Dirac formalisms, are introduced and compared. The various approaches, numerical and analytical, to calculating the coupling within these formalisms are discussed. Finally, the merits and disadvantages of these methods, and those of more complete ab initio treatments, are given.

Calculating Exchange Interactions
The IEC between a pair of magnetic objects is usually defined as the energy difference between the parallel and antiparallel alignments of the localised magnetic moments on the objects. We distinguish the indirect interaction by assuming that it is mediated entirely by the conduction electrons of the host material, in contrast to the direct exchange which depends on the overlap of the impurity electron orbitals and which decays very abruptly. To begin, we consider an initial configuration of two magnetic impurities at sites A and B, separated by a distance D, in the host medium. We introduce a Green function,ĝ, that describes the electronic structure of this system. For the moment we will keep the form of this Green function, and its associated Hamiltonian, as general as possible. We assume that it describes in sufficient detail electron propagation in the host material and also spin-dependent potentials localised at the magnetic impurities, which reflect the fact that electrons in the system are subject to a different potential in the magnetic regions than elsewhere, as shown schematically in Figure 2. The resultant potential profiles seen by "up" and "down" spin electrons are the result of solving the Hubbard Hamiltonian and are illustrated by red and blue lines respectively. The allowed energies of the system are quantised and furthermore, changing the separation between the magnetic objects shifts the allowed energies relative to the Fermi energy of the system giving rise to an oscillatory coupling [52]. The coupling can be calculated by summing over all the energy levels below the Fermi energy and taking the difference between the cases with spin potentials corresponding to moments aligned parallel or anti-parallel. In this manner, the calculation of the IEC reduces to the calculation of an energy difference between two distinct configurations. Such a calculation can be simplified using the Lloyd formula [51], which allows the calculation of the energy difference directly without the need to calculate the total energy of either configuration.
The Lloyd formula expression to calculate the total change in the energy of a system due to a perturbation is given by whereĝ as above is the Green function describing the unperturbed system,V is the applied perturbation potential and f (E) is the Fermi function. The unperturbed system consists of two moments embedded at sites A and B, and aligned parallel along the z-direction so that the angle between them, θ = 0. This setup is shown schematically in the panel a of Figure 2. For simplicity we consider one magnetic orbital on the impurity sites, but the calculation can be easily generalised to a more realistic magnetic impurity. We associate a bandcentre (δ) and exchange splitting (V ex ) with this orbital as illustrated in Figure 2. We now introduce a spin perturbation which rotates the magnetic moment at B by an angle θ with respect to that at A. This perturbation is given by whereσ z andσ x are the relevant Pauli matrices. Since the magnetic moments break the spin degeneracy of the electrons, the Green functions must be written in terms of their up-and down-spin components. In the initial collinear configuration there is no mixing between the spin bands and g BB is diagonal in spin space. It can be shown that where g ↑ AB (g ↓ AB ) are the off-diagonal elements of the up (down) spin Green function connecting sites A and B. Taking θ = π, corresponding to an antiparallel alignment of the moments we find, for zero temperature, that where we have chosen the sign convention where a negative value of the coupling, J, corresponds to a preferential parallel alignment of the moments and a positive value to a preferential antiparallel alignment. When performing the integral in Equation (4) numerically we can rewrite the integral over the imaginary axis where the integrand tends to be smoother and easier to integrate. With an imaginary axis integration, the expression for the coupling becomes We note that in these expressions the distance dependence is entirely contained within the product of off-diagonal Greens function elements and the behaviour of these quantities will dictate the behaviour of the interaction as the separation between the moments is varied. The expression for the coupling given in Equation (4), with the log term in the integrand, is exact. That is, it gives the exact energy difference between the FM and AFM configurations of the system described by the Green functions used. In the next section we will discuss how the relevant Green functions will be calculated, but first we focus on a common perturbative approach to calculating the coupling. To proceed analytically it is useful to note that Equation (4) can be written as a perturbation expansion in powers of the exchange splitting V ex and when expressed to leading order in V ex gives where g BA (E) is the spin-independent Green Function describing electron propagation in the pristine host material [53]. This expression for the coupling is equivalent to the Ruderman-Kittel-Kasuya-Yosida (RKKY) approach, initially developed to describe the coupling mechanism of nuclear magnetic moments [40], then expanded to describe a wider range of indirect coupling phenomena [41,42] and generalised to provide a model capable of reproducing experimental observations [43,44]. This treats the coupling as a consequence of the spin polarisation of the conduction electrons of the host by the magnetic objects. Under this approach, the coupling is written as an effective direct coupling, J BA s B · s A , between the two moments, where the coupling strength is given by where λ is an adjustable parameter representing the magnitude of the coupling between localised spins and conduction electrons and χ is the static magnetic susceptibility which relates the response of the host magnetisation to a static magnetic field. Writing this quantity in terms of Green functions, and the equivalence of Equations (6) and (7) becomes clear. Thus the RKKY approach can be considered a second-order perturbational approximation to the IEC. The approaches are in general equivalent, except when V ex >> 1 or when the host material possesses certain symmetries [53][54][55][56]. The RKKY approximation generally provides a good description at large separations when the coupling is quite small. An important contrast between the Lloyd formula approach and the RKKY approach is the choice of Green function used. The former approach avails of the spin dependent Green functions in the ferromagnetic configuration, g ↑ BA and g ↓ AB whereas the latter approach generally uses their pristine, spin-independent counterparts. We note that use of the pristine GFs essentially sets the value of the bandcentre, δ, to zero. An intermediate approach can be taken by using the spin-dependent Green functions with an RKKY-like expansion of the logarithm in the Lloyd expression. The general form expected for the coupling is where k F is the Fermi wavevector in the conducting host and the exponent α determining the rate of decay with separation depends on the dimensionality of the system. Most of the analytic investigations into the IEC in graphene make use of the RKKY approximation. Although we have written it here in terms of energy dependent Green functions, it is often formulated differently. A common alternative formulation is in terms of time dependent Green functions, the Fourier transforms of the GFs used here, although conflicting results have been reported depending on whether the real or imaginary time representations are used [57][58][59]. Numerical approaches to indirect exchange interaction generally proceed in one of two ways. The first uses numerical methods to calculate either perturbative or exact expressions for the coupling like those given above. Another method involves direct calculation of the energies of the individual FM and AFM configurations. This is a common approach in ab initio calculations [60][61][62][63][64] although it has been employed also for tight-binding models [65,66]. The computational cost of this type of calculation tends to limit these investigations to small or medium values of separation and it can be difficult to perform a meaningful comparison with results for the asymptotic behaviour calculated in other ways.

Calculations in Graphene
Now that we have introduced the general methodology to calculate indirect exchange couplings, we shift our focus to examine the specific case of interactions in graphene. The electronic structure of graphene plays a key role in determining the interaction properties and so we first outline the different methods used to calculate or approximate the band structure. We then discuss how the Green functions required for some coupling expressions can be found and compare the different approaches used to calculate the coupling.
The electronic structure of graphene near the Fermi energy is principally determined by π bonds formed from the out of plane 2p z orbitals of carbon atoms situated on the vertices of a hexagonal, or honeycomb, lattice. The remaining carbon orbitals hybridize in plane to form very stable σ bonds which give graphene its extraordinary strength. Thus a single-orbital nearest-neighbour tight-binding Hamiltonian is usually more than adequate when describing the electronic properties of graphene. To write this Hamiltonian we make use of the two-atom unit cell shown in Figure 1. Thus every site in the lattice is denoted by a vector r giving the unit cell location and an index n = 0, 1 differentiating between the two sites within the unit cell, illustrated by filled and hollow symbols in the figure. The tight-binding where |rn⟩ represents the orbital at the site given by {r, n} and t n,n ′ r,r ′ is the electronic hopping term between two such orbitals. The sum is restricted to orbitals at neighbouring sites and the non-zero hopping terms all take the value t = −2.7 eV. One consequence of only allowing nearest neighbour hoppings is that graphene is a bipartite lattice. The sets of filled and hollow symbols in Figure 1 can be regarded as intersecting triangular sublattices that together form the honeycomb graphene lattice. Introducing hopping terms beyond the first nearest neighbour break the bipartiteness of graphene. However, since these terms are much smaller than the nearest neighbour hoppings, graphene behaves effectively as a bipartite lattice. It is this feature of the graphene lattice, combined with the half filling of its π orbitals, that determine many of the features of the RKKY interaction [57].
Using Bloch theorem methods to take advantage of the lattice periodicity, the associated eigenvalues and eigenvectors of Equation (10) are found to be and 3 . The bands given by Equation (11) are shown plotted along the high symmetry points of the graphene Brillouin zone by the solid lines in Figure 3a.
We note that since the p z orbitals in carbon contain one electron each, then the band is half full and the Fermi energy of the undoped system is E F = 0. This is exactly the point where the two bands given by ϵ − and ϵ + touch and the DOS vanishes. Graphene is thus a zero-bandgap semiconductor, or a semi-metal. The resultant Fermi surface consists of six discrete points, only two of which are unique, lying at the vertices, or K points, of the hexagonal first Brillouin Zone (BZ) illustrated in Figure 3. Many of the interesting properties belonging to graphene arise due to the shape of the bands near the Fermi energy. Unlike in conventional semiconductors where the bands are parabolic, the band structure of graphene is linear near E F . This results in electrons or holes near the Dirac points having zero effective mass and behaving like relativistic particles which can be described using the Dirac equation from Quantum Electrodynamics. It is worth investigating briefly how graphene electrons are approximated in the linear regime and the range of energy values over which the approximation is valid.
The region of the graphene spectrum around E F = 0 can be approximated linearly to simplify the calculation of physical properties. For a point in reciprocal space near the Dirac point, k = K + δk, we find where the Fermi velocity of graphene, v F = 1 dE dk ∼ 10 6 m s −1 . In Figure 3a we compare the linear approximation of the band structure (dashed lines) with that calculated previously (solid lines). It is seen to be a very good approximation in the region surrounding E = 0, but loses accuracy quickly outside this regime.
To calculate the coupling using the expressions in the previous section, we need the Green function connecting two sites on the pristine graphene lattice. There are a number of different ways to calculate this quantity and use it to calculate the coupling, but for illustration we will focus on a real space, energy dependent Green function calculated from the dispersion relation and eigenvectors of the non-linearised tight-binding Hamiltonian. The Green function we require is given by an expression of the form where and the reciprocal space integral is over the first BZ, where for convenience, we can define a rectangular BZ consisting of segments drawn from multiple neighbouring BZs whose area equals that of the hexagonal BZ. A possible rectangular BZ is shown by the shaded area in Figure 3b. We note that the Green function expression takes a different form for sites on the same or on opposite sublattices. The most obvious, if not most efficient, way to proceed from here is a brute-force two-dimensional numerical integration. However, we have shown previously [67] that it is possible to perform one of the integrals completely analytically, greatly reducing the numerical cost of the Green function calculation. Furthermore, for separations greater than a few lattice spacings in the high symmetry zigzag and armchair directions it is possible to approximate the remaining integration to a high degree of accuracy using the Stationary Phase Approximation (SPA), to yield a closed form expression for the Green function where A(E) is an energy-dependent coefficient and Q(E) can be identified with the Fermi wavevector in the direction of separation. The full analytical expressions for these quantities for armchair and zigzag separations can be found in [67]. Once calculated using any of these methods, the Green functions can be substituted directly into the RKKY integration in Equation (6). Alternatively, they can be used to calculate the spin-dependent Green functions required for the non-perturbative total energy expression in Equation (4), where the Dyson equation is used to add a spin-dependent potential at the magnetic sites. At this stage, it is worth mentioning some of the alternative methods used to calculate the graphene Green function for RKKY calculations in the literature. Many of the earlier, analytic studies employ the linear approximation to the graphene band structure [57,58,[68][69][70] and may also use time-dependent Green functions in lieu of the energy dependent form shown here [57][58][59]. There is some discrepancy between early results for both the sign and decay rate of the interaction (e.g., [57,58,69]), but many of the key features to be discussed in the next section are present in [57] and later confirmed by other studies. Numerical investigations using Green function techniques have tended to use an energy dependent form, but calculated in a variety of ways [67,71,72]. Sherafati and Satpathy [71] employ both full numerical integration of the Green function integral in Equation (14) and the recursive Horiguchi method, which relates the graphene Green function to those of the triangular lattice, which in turn are expressed in terms of elliptic integrals [73]. Although more efficient than the full integration, the Horiguchi method was only stable for small values of separation [71]. In the same work, an analytic treatment in the linear dispersion regime using energy dependent GFs confirmed, with some refinement, many of the features predicted previously using time dependent GFs. In a study of RKKY interactions in disordered graphene sheets, Lee et al. [72] used the Kernel Polynomial Method [74] to express the Green functions in terms of Chebyshev polynomials with recursively calculated coefficients. Working outside the Green function formalism the interaction can be calculated by determining the total energies of the FM and AFM alignments of the moments. Such an approach can use tight-binding models similar to those presented here [65,66] or more complete ab initio calculations [60][61][62][63][64]. Such methods often rely on periodic boundary conditions and relatively small unit cells, leading to problems like interference between moments in neighbouring cells [75] and difficulty in calculating the interaction for separations greater than a few multiples of the lattice parameter.

Interaction between Simple Magnetic Impurities
In this section we consider localised magnetic moments that are associated with a single lattice site in graphene. These can be thought of as substitutional magnetic impurities replacing a single carbon atom in the lattice, or as moments induced at a lattice site by the presence of an adsorbate. This type of impurity is the most commonly considered in the literature and the exchange interactions between pairs of these impurities have been investigated using many of the methods discussed in the previous sections. The expressions for the couplings given in Equations (4) and (6) are easily applied to moments of this kind as the Green functions required are either those for pristine graphene with a simple spin-dependent perturbation at the magnetic sites, or without this perturbation if the RKKY expression is used. Throughout this section we will consider two magnetic impurities occupying sites A and B in the graphene lattice, separated by a distance D, with an indirect exchange interaction J AB .
In the earliest studies, the RKKY interaction in graphene was often considered in conjunction with Friedel Oscillations (FO). The latter involve distance dependent charge fluctuations relative to an impurity site and there is a large degree of overlap in the physics involved. Using an analytic Dirac linearisation approach to study FO, Cheianov and Fal ′ ko predicted that charge density FO would decay as δρ ∼ D −3 , but the RKKY interaction as J ∼ D −2 [69]. Wunsch et al., meanwhile, found the magnitude of induced moments to decay as D −3 away from a magnetic impurity and predicted the RKKY interaction to behave likewise [76]. Before these reports, an early study by Vozmediano et al. investigating possible ferromagnetism in graphene layers arising from localised moments near lattice distortions had examined the RKKY interaction between moments surrounding elongated cracks in graphene [68]. Although more complicated structures than considered elsewhere, a decay rate of D −3 for a non-oscillatory, always ferromagnetic interaction was reported. An important paper by Saremi clarified many of the features of the interaction [57]. Firstly, a general proof was given regarding the sign of the RKKY interaction in a half-filled bipartite lattice. This showed that, due to electron-hole symmetry in such systems, the coupling was ferromagnetic between moments on the same sublattice and antiferromagnetic between moments on opposite sublattices. In the same work, analytic calculations were performed in the linear dispersion regime to determine the functional form of the RKKY interaction. These calculations included the explicit use of a cut-off function to prevent divergences arising from high energy contributions. There has since been some debate about the form, or even the necessity, of the cut-off function and of the nature in which the time-dependent Green function integrations are performed [57][58][59]65,67,71,77]. However, many of the principal features of the interaction reported by Saremi have been verified by other studies, including a contemporaneous paper by Brey et al. [70]. Bunder and Lin suggested that differences between the real time and imaginary time approaches could lead to a breaking of the theorem regarding the sign of the interaction [58]. In particular they report sign changing oscillations for the interaction in zigzag-edged graphene nanoribbons. However, fully-numerical calculations of the coupling by Black-Schaffer do not find such sign changes, which the author attributes to the perturbative nature of the previous calculations [65]. In the same work, the main predictions made by Saremi are verified using a numerical exact diagonalisation method to find the energy difference between the FM and AFM orientations. A minor refinement is the introduction of a π phase shift in the cosine term for moments on opposite sublattices. In fact, Sherafati and Satpathy showed that this phase shift is not a constant but depends on the direction of separation, taking a maximum value of π for zigzag separations and vanishing for armchair separations [71]. This additional phase shift term can be added to the three features proposed by Saremi to provide a complete description of the RKKY interaction between simple single-site magnetic moments in graphene.
Using the stationary phase form of the Green function in Equation (15), we have previously demonstrated how the principal features of the interaction in the RKKY approximation can be calculated [67]. The integration over energy in Equation (4) can be reduced to a sum over Matsubara frequencies. The functions B(E, ε) = A 2 (E, ε) and Q(E, ε) are expanded around E F and in the low temperature limit T → 0, we find where is the distance-independent coefficient for the ℓ-th term in the power series, ℓ is a non-negative integer and B (ℓ) (E F ) is the ℓ-th order energy derivative of B(E) evaluated at E F , resulting from its Taylor expansion. In general the leading term in the series should determine the asymptotic decay rate of the coupling. For the undoped case it can be shown that the coefficient B (0) (0) = 0, so that the l = 1 term dominates and J(E F = 0) ∼ D −3 . The exact oscillatory term predicted by Saremi is reached by forcing the real part of the Green function to vanish at E F = 0. Figure 4 shows a numerical calculation of the coupling, calculated using Equation (5), for substitutional magnetic impurities in graphene. The coupling is plotted as a function of separation for both armchair and zigzag directions, and for impurities on the same and on opposite sublattices. Results are shown for both undoped (E F = 0.0) and doped (E F = 0.1|t|) cases. The full range of features discussed in this section are present. In particular we note the D −3 decay rate and absence of sign-changing oscillations in the undoped cases. We note the presence of a period-3 oscillation in the zigzag direction compared to the monotonic decay for the armchair case. This results from slightly different commensurability effects in the two directions. In undoped graphene, the oscillations are determined by cos 2 (k F · D). The components of the Fermi wavevector in the armchair and zigzag directions are 2π √ 3a and 2π 3a respectively. Since the spacing between unit cells is √ 3a for armchair and a for zigzag separations, the cosine squared term takes a constant value of 1 in the armchair direction but cycles repeatedly through the sequence 1 4 , 1 4 , 1 as the separation is increased in the zigzag direction. Also evident in the zigzag separation plot is the π phase shift between same and opposite sublattice cases. When the system is doped, as in the bottom panels, a longer-ranged interaction decaying as J ∼ D −2 is seen. Long-wavelength, sign-changing oscillations are also introduced by the change in the shape of the doped Fermi surface.  Armchair separation Zigzag separation

Beyond the Simple Case
In the previous section we highlighted the features of the magnetic interaction in graphene for the simplest kind of impurities. These impurities are the kind usually investigated within the RKKY approximation. In this section we consider a number of extensions to the simple Hamiltonian describing the system and examine the effects they have on the coupling. We consider first a more general or realistic description of the magnetic impurities and examine the effect of impurity parameterisation on the interaction. We show that the use of the pristine Green function in the RKKY approximation ignores these effects which can in principle overturn some of the results noted in the previous section. Then we move to the case when the impurity connects to more than one atom in the graphene lattice and examine the effect on the decay rate of the interaction. Finally, more detailed descriptions of the graphene electronic structure are considered. We discuss the possible effects of electron-electron interactions and spin-orbit coupling and examine the results of some ab initio calculations.

Impurity Parameterisation
Calculations using the RKKY approximation do not generally account properly for the local spin dependent potentials describing the magnetic moments. In fact, the parameterisation of the magnetic impurity is generally neglected through the usage of the pristine Green functions. In this section we consider the effect of three parameters that can be used to characterise the magnetic impurities. In addition to the magnetic moment (m) and band-centre shift (δ), shown in Figure 2 and which can be calculated, for example, using self-consistent mean-field calculations, we also consider the hopping potential between the lattice carbon sites and the impurity site (t ′ ) which should differ from the carbon-carbon hopping. This approach is similar to the Anderson model describing localised magnetic impurity states in metals [78]. These parameters will vary between the different magnetic species that can be chosen as the embedded impurities.
In Figure 5 we plot the coupling between substitutional impurities on the same sublattice, calculated using the Equation (4) and exact Green functions, as a function of separation along the armchair direction for three different parameter sets (m, δ, t ′ ). The first of these (0.6, 0.0, t) closely replicates the results of the RKKY approach as it considers only a band-splitting, has no band-centre shift and uses the carbon-carbon hopping value. The numerical calculations in the previous section used a similar parameterisation. The middle and bottom plots use the parameters (0.6, 5.0, 0.8t) and (0.6, 8.0, 0.6t) respectively. In these plots we see the formation of an unusual feature not predicted by the RKKY approximation. For quite a large range of distances we note a preferential antiferromagnetic alignment between the moments before the sign flips and the standard ferromagnetic coupling with a 1/D 3 decay is recovered. A similar sign-changing behaviour has been reported in nanotubes [49] and also by ab initio calculations attempting to probe the interaction in graphene [64]. It is worth examining further how this feature depends on the parameterisation of the magnetic moments. In Figure 6 we present a number of phase diagrams showing the sign and strength of the coupling for different values of these parameters. Each diagram represents an area of (m, δ) phase space with ferromagnetic (antiferromagnetic) couplings given by a blue (red) shading that is darker for larger magnitude couplings. The diagrams on the top row correspond to a separation of 10 √ a between the magnetic moments and from left to right show the cases of t ′ = 1.0t, 0.8t, 0.4t. The bottom panels show the same cases for a larger separation of 40 √ a. By examining the border between the blue and red regions in these plots we can infer under what circumstances the sign-change behaviour described above occurs. In all cases the border position varies only weakly with the magnetic moment (m) or hence the band-splitting (V ex ). A stronger dependence is found on the band-centre shift (δ) and we find that, in general, an anti-ferromagnetic alignment is found above a critical value of δ. The band-centre shift is strongly dependent on the occupation of the magnetic orbital. For a bipartite lattice like graphene a substitutional impurity with a half-filled orbital gives zero band-centre shift when t ′ = t. As we move away from half-filling a larger band-centre shift is required to return the correct band occupation. For smaller values of t ′ we note that the border between the blue and red regions shifts towards the left, meaning that smaller band-centre shifts will lead to an AFM alignment. Increasing the distance between the impurities reduces the phase-space area corresponding to an AFM alignment by shifting the border to the right and requiring larger band-centre shifts. This finding agrees with the distance-dependent plots in Figure 5 for fixed parameters which show that in the asymptotic limit the coupling changes sign and returns the FM alignment predicted by the RKKY approximation. However, as in the bottom panel of Figure 5, the magnitude of the coupling has essentially decayed to zero before the sign change occurs so the only significant coupling between two such impurities is antiferromagnetic. At this point, we note that the theorem proposed by Saremi [57] suggesting that the same (opposite) sublattice coupling is always FM (AFM) only applies for the case of bipartite lattices with half-filling. The configurations discussed here break the half-filling requirement and electron-hole symmetry, at least locally, and this has a dramatic effect on the nature of the interaction. In particular, depending on the moment parameterisation, a strong antiparallel alignment between two substitutional magnetic impurities in graphene may persist to considerable separations. This result is contrary to many predictions made regarding the RKKY interaction in graphene and shows that caution needs to be applied when extrapolating the results of simple calculations to more realistic magnetic impurities.

Multi-Site Impurities
The features of the RKKY interaction discussed in Section 3 contain a strong sublattice dependence. A substitutional or top-adsorbed impurity is associated with a single site on the graphene lattice and therefore also with a single sublattice, as illustrated in Figure 7a. We have seen how the sign of the interaction, at least for simple impurity parameterisations, depends on whether the interaction is between sites on the same or on opposite sublattices. However, another type of impurity can be considered which is associated with an equal number of sites on both sublattices. Two examples of this type of impurity are shown in Figure 7. These are the bridge impurity (b), adsorbed onto the lattice midway between an atom from each sublattice, and the plaquette, or centre-adsorbed, impurity at the centre of hexagon (c) and bonding equally to the six surrounding carbon atoms. Plaquette impurities have been investigated more thoroughly than bridge impurities and a number of interesting predictions have been made. An important difference is noted between the results for incoherent and coherent type moments [57]. For incoherent moments, the calculation essentially consists of adding the contributions of the thirty-six single-site interactions between the associated sites of one moment and the other. Coherent moment calculations are more complete and can be performed, for example, using Equation (4) with the spin-dependent Green functions of a graphene sheet with the two impurities connected to the graphene sheet using the Dyson equation and an appropriate perturbation potential. Saremi noted a decay rate of D −3 for incoherent plaquette impurities [57], a result confirmed by further studies [65,70,71]. However, Saremi also noted that for the coherent case the D −3 term vanishes and that the decay in this instance is faster. This is similar to the decrease in the interaction range reported for plaquette impurities in metallic carbon nanotubes [48], where the decay rate was found to be D −5 compared to the D −1 rate predicted for substitutional or top-adsorbed impurities. Uchoa et al. have since reported a decay rate of D −7 for plaquette impurities in graphene [79].
The sign of the interaction between plaquette impurities is also of interest. An always AFM coupling is widely reported in the literature [57,65,71], although Uchoa et al. report a FM coupling for smaller values of separation [79]. The prevalence of an AFM coupling can be understood from the single-site coupling results. We may initially expect the plaquette impurity to average out the sublattice-driven preference for an always FM or AFM coupling. However, we recall that the magnitude of the AFM coupling in the single site case is three times larger than that of the FM coupling. Thus the AFM contributions from inter-sublattice terms are larger than FM contributions from intra-sublattice terms and an overall AFM coupling is found between plaquette impurities. Another interesting feature for this impurity type is that the period-3 oscillations noted for zigzag separations are no longer present and a smooth monotonic decay of the coupling is noted for separations in all directions. The jagged features are averaged out by the fact that a plaquette impurity is associated with an equal number of sites from each group defined by the period-3 sequence noted for single-site impurities. Sherafati and Satpathy [71] have examined the case of bridge impurities and find a rare example of sign-changing oscillations in the interaction in undoped graphene. They report that every third value of separation has an AFM coupling, and the other two are FM.

Spin-Orbit Effects, Electron-Electron Interactions and ab initio Calculations
Having spent some time examining the effects of more detailed impurity parameterisations, it is now worth considering the role a more complicated description of the graphene electronic structure may play. Such a description may include additional terms in a model Hamiltonian calculation or a full ab initio treatment using a Density Functional Theory (DFT) approach. An early study of this type was undertaken by Dugaev et al. [80] and investigated the interaction in a graphene system where a small gap was introduced by the inclusion of spin-orbit effects. Using an analytic approach based on the Dirac equation, the prediction of an always FM interaction with a long tail exponential decay was made.
Black-Schaffer [66] included the effect of electron-electron interactions in a numerical study of the RKKY interaction. These interactions were included using a mean-field solution of a Hubbard Hamiltonian with a finite onsite Coulomb repulsion term U at each site in the graphene lattice. Similar methods have been used previously to investigate intrinsic magnetism in graphene arising from the spin polarisation of localised states near zigzag edges. A number of interesting effects are noticed as the value of U is increased. The oscillatory features in the zigzag direction are removed, as is the factor of three amplitude difference between inter-and intra-sublattice couplings. Furthermore, the decay exponent decreases from 3 for U t = 0 to below 2 when U t > 2. At this value of U the interaction is also isotropic, with no difference between the interactions for armchair and zigzag separations.
A few ab initio studies have been performed which address exchange couplings between localised magnetic moments in graphene. These studies are generally difficult to compare with the other results presented here for a number of reasons. Due to a much larger computational cost, ab initio calculations are generally only performed for very small values of separation. The majority of these calculations use periodic boundary conditions, which can result in magnetic interactions between impurities in neighbouring unit cells [75,81]. In fact, such interactions can lead to unusual results in ab initio studies examining a single magnetic impurity due to the non-oscillatory nature of the RKKY interaction in graphene. Preferential AFM alignment between impurities in neighbouring unit cells can potentially lead to the total suppression of the magnetic moment [31,75]. Periodic unit cell calculations within a tight-binding scheme can minimise this effect by including a large "padding" region around the impurities [65], but this is usually not feasible in DFT. Signatures of the RKKY interaction can also be seen in DFT calculations when the effective distance between moments in single impurity concentrations is varied by increasing the size of the unit cell [32,60]. From a study of pairs of localised moments arising from simple vacancy defects, Pisani et al. [61] extracted a decay rate of D −1.43 for moment separations of up to 25Å. This agrees with Black-Schaffer's calculations including electron-electron interactions when the onsite Coulomb repulsion term is U t = 2.1. A study by Santos et al. [62] however extracts a decay rate of D −2.43 for a FM decay between substitutional cobalt impurities on the same sublattice. This result is significantly closer to the D −3 rate predicted without electron-electron interactions. Calculations between Fe atoms meanwhile suggest a preferred AFM alignment between impurities on the same sublattice [75], suggesting that the nature of the interaction depends strongly on the magnetic species chosen, in agreement with the findings of Section 4 of the current work. Moments arising in hydrogenated armchair-edged graphene nanoribbons generally preferred an AFM alignment [63]. Carbon adatoms, adsorbed in the bridge configuration, were found for small separations to have a FM alignment. For larger separations both FM and AFM alignments were found for different values of separation [64].
These findings show that there is still quite a large amount of discrepancy between the interactions predicted using simple RKKY-like models and more complete DFT calculations, which appear very strongly dependent on impurity configuration. Furthermore, direct comparison of results from both methods is often thwarted by the separation constraints imposed by DFT calculations and the possibly oversimplistic treatment of magnetic impurities in model calculations.

Manipulating the Interaction
In this section we discuss a number of factors which may allow manipulation of magnetic interactions in graphene. A degree of control over the properties of the interaction may be required for spintronic applications, for example to switch on and off spin currents or to dynamically change the magnetic ordering in a system. Another important consideration is that the interactions discussed so far have been very short-ranged, due to a fast decay rate arising from the peculiar electronic dispersion relation in graphene. Methods of augmenting the interaction, or slowing the rate of decay, may thus help overcome some of the difficulties in realising magnetically-doped graphene devices.

Geometry Effects
A number of studies have looked at magnetic interactions in graphene systems with geometries other than the standard two-dimensional sheets considered so far in this work. These include nanotubes [47][48][49][50]82,83], nanoribbons [58,63,65,66,84,85], nanoflakes [86] and multilayer graphene systems [59,[87][88][89][90]. Quasi-one-dimensional graphene systems, such as nanotubes and nanoribbons, are expected to display longer ranged interactions than those present in two-dimensional graphene sheets. Indeed, metallic nanotubes have been shown to display a long ranged D −1 interaction when substitutional, top-adsorbed or bridge-adsorbed impurities are considered [47,48]. However, similar to the case of graphene, the interaction between plaquette impurities is found to decay at the much faster rate of D −5 . The sublattice dependent sign rules discussed for graphene also hold for nanotubes within the RKKY approximation, but can be broken when the parameterisation of the impurities is altered [49]. The eigenstates of carbon nanotubes are subject to periodic boundary conditions in the circumferential direction which ensure the equivalence of all lattice sites. This is not the case in nanoribbon systems where edges are present and each lattice site can be characterised by its distance from the edges. Thus we should expect the interaction to contain a dependence on the location of the impurities across the width of the ribbon. Such a dependence is found in many properties of doped ribbons, including binding energies [91], conductance [92,93] and the magnitude of magnetic moments [94]. The presence of localised edge states for zigzag edge geometries [17][18][19] should also be expected to introduce features not present in graphene sheets or nanotubes. Black-Schaffer [65] examined the interaction in nanoribbons, finding that armchair edges did not significantly alter the coupling from the bulk case. However, for localised moments at the edge of zigzag ribbons, an exponentially decaying interaction was reported. The interaction was FM for moments along one edge and AFM for opposite edges. However, the magnitudes of the FM and AFM interactions were equal, unlike the larger AFM magnitude found for the bulk case. Towards the ribbon centre, the results for the bulk were recovered. Including electron-electron interactions in the calculation suggested that the distance dependence vanished for reasonable values of U , leading to an extremely long-ranged interaction [66]. Szałowski studied the RKKY interaction in finite graphene nanoflakes, finding that for certain geometries the addition of a single additional charge carrier could change the sign of the interaction [86]. A recent work by Klinovaja and Loss extends the study of the interaction in one-dimensional graphenes to include spin-orbit interactions [95]. The RKKY interaction has been studied in bilayer graphene systems by several authors [59,[87][88][89]. The quadratic nature of the dispersion relation in bilayers means that many of the features present in monolayer graphene are removed. Another interesting feature is that, for A-B stacking, a plaquette impurity in one layer is embedded above a single graphene lattice site on the other layer. Extending the study to more than two layers, Jiang et al. [90] find a number of different decay rate possibilities, depending on the positions of the two impurities and on the number of layers.

Doping and Disorder
The key features of the RKKY interaction when the Fermi energy is shifted away from half-filling are clear from the bottom panels of Figure 4. These are namely a longer ranged interaction that decays as D −2 and sign-changing oscillations, the periods of which depend on the component of the Fermi wavevector in the direction of separation. The appearance of oscillations is due to the breaking of the commensurability effect between the Fermi wavevector and the lattice spacings. Many of the studies of the IEC in graphene establish these key features in doped or gated graphene, but they have been explored in more detail by Sherafati and Satpathy [96], where the possibility of controlling the sign of the interaction with a gate voltage is raised. The possibility of controlling the interaction using a gate voltage has also been discussed for bilayer graphene [88].
Lee et al. [72] examine the effects of nonmagnetic disorder on the RKKY interaction in monolayer graphene. The Anderson model was employed to study both onsite and hopping term disorder. Onsite disorder was found to induce sign-changing oscillations in the coupling in certain directions, whereas off-diagonal disorder only effected the amplitude, and not the sign, of the interaction. This suggests that sign-changing behaviour is introduced by the breaking of sublattice symmetry, in a similar manner to that discussed for impurity parameterisation in Section 4. For weak disorder, the interaction retained similar decay behaviour to the clean case. In the strongly disordered, localised regime exponential suppression of the coupling with increasing moment separation was observed. In a recent work, the same authors examine the interplay of disorder and gate voltage [97].

Strained Graphene Systems
Strain-engineering of graphene systems has received a lot of attention in recent literature due to the possibility of tuning many of the physical properties of graphene [98][99][100][101][102][103][104][105][106][107][108][109]. The degree of tuning is enhanced by the different types of strain that can be applied. Apart from simple uniaxial strains [98,110], more exotic features like creases and bubbles can be introduced [101,104,106,[111][112][113]. In a recent work, we have explored the possibility of manipulating the IEC in graphene by the application of uniaxial strains parallel or perpendicular to the moment separation direction [109]. Another recent work also explores this topic [108]. We find that both amplification and suppression of the magnetic coupling can be achieved. We reported a general trend of amplification for strain perpendicular to the moment separation direction, or suppression for strain parallel to this direction. Also noted are oscillations in the amplification as strain is increased for moments separated in non-armchair directions. Such oscillations suggest the intriguing possibility of selectively turning on or off the coupling between moments and also of controlling the inter-and intra-sublattice couplings independently. Since the IEC underpins physical features, including overall moment formation and magnetotransport response, the ability to fine tune the coupling with strain may lead to interesting spintronic applications. Thorough investigation of graphene systems with magnetic impurities and with different types of strain applied may yield a further range of tuneable spintronic properties.

Spin Precession and Dynamic RKKY
One of the major obstacles in the path of exploiting the RKKY interaction between magnetic impurities in graphene for spintronic devices is the very short range of the interaction. It has been suggested that the range of RKKY-like interactions can be augmented if the magnetic moments are set in motion, for example, if they are set to precess by the application of a suitable time-dependent magnetic field [114,115]. Such an interaction is driven by non-equilibrium spin currents emanating from the precessing moments [116]. The magnitude of the interaction can be measured by a quantity called the dynamic spin susceptibility, which describes the response of the magnetism of the system to a dynamic magnetic perturbation. Investigation of this quantity in carbon nanotubes [82] has revealed a decay rate slower than that of the static coupling. Further studies of spin dynamics in graphene systems have suggested the use of these materials as spin waveguides [117], spin-pumping transistors [118] and spin current lenses [119]. Recent experimental evidence also suggests possible long-range spin current behaviour in graphene [120]. In a recent study [121], we examined the spin susceptibility in graphene as a dynamic analogue of the static RKKY coupling, with a particular focus on the separation dependence of the interaction. We find that the decay rate for the dynamic RKKY approaches D −1 for both doped and undoped systems, giving a significant augmentation of the interaction range. Furthermore, the dynamic interaction can be related to fluctuations in the lifetimes of magnetic excitations at the impurity sites. Such quantities have been probed in other materials using inelastic scanning tunnelling spectroscopy [122][123][124][125], suggesting that this method may be used to detect signatures of the dynamic RKKY interaction in graphene.

Conclusion and Experimental Considerations
In this review we have examined many aspects of indirect exchange interactions in graphene systems. The various approaches to calculating the coupling, ranging from analytical solutions of the RKKY approximation using the Dirac approximation, to numerically expensive ab initio total energy calculations, via tight-binding energy difference calculations with varying levels of impurity parameterisation, were discussed. Each of these methodologies has merits and pitfalls, and the contrasting results suggest that there is still some work to be done before useful comparison can be made with experimental results. On the experimental side, it is understandable that with a decay rate as fast as D −3 it is difficult to probe the interaction for any reasonable separation. The presence of magnetism in disordered graphene systems may indicate the presence of an exchange coupling between magnetic moments formed around defects. Nuclear magnetic resonance experiments reveal that these defects have indeed magnetic moments, since they couple to implanted Fe atoms [126]. However, whether or not these moments couple with each other, or with the graphene lattice, to form a ferromagnetic state is a controversial subject and many of the results in this area have proved difficult to reproduce [5]. Progress in spin-dependent scanning tunnelling spectroscopy has enabled the RKKY interaction to be measured directly from the magnetisation curves of individual magnetic adatoms [127,128]. As discussed above, signatures of the dynamic analogue of the RKKY interaction may be detectable from the magnetic excitation lifetimes of individual moments measured using inelastic scanning tunnelling spectroscopy [122][123][124][125]. In this case, the longer range of the interaction may make experimental detection more feasible [121]. Using a gate voltage to increase the interaction range may also be necessary to probe the interaction in the static case. To our knowledge these types of experiment have not yet been performed with a graphene host medium. Signatures of indirect exchange interactions may also be present in transport measurements taken in magnetically-doped graphene systems [49,129], since different orientations of the magnetic moments lead to different scattering regimes for up-and down-spin electrons.
Future work on the topic of indirect exchange interactions in magnetically doped graphene should focus on the gaps and discrepancies between the various theoretical models, and between these models and possible experimental configurations. In particular, the strength of the interaction for specific magnetic impurities and whether it survives from the ideal model to experimentally realisable conditions of finite temperature and lattice disorder are important considerations. Another area of exploration is the interplay between magnetic impurities and other features that can be introduced into the graphene lattice. The topics of different geometries, in the form of nanoribbons and nanotubes, and of geometrical distortions, in the form of strain, were briefly discussed in this review and results were mentioned for simple cases. More complicated geometries or distortions may introduce new features and more possibilities to manipulate the interaction. In moving towards technological application of RKKY-like interactions, the study of larger ensembles of magnetic impurities will be necessary. The magnetic and electronic properties that emerge from such ensembles should be dictated by the pairwise indirect exchange interactions discussed in this review, and may lead to a range of interesting features and applications.