# The RKKY Coupling between Impurity Spins in Graphene Nanoflakes

###### Abstract

We calculate the indirect charge carrier mediated Ruderman-Kittel-Kasuya-Yosida (RKKY) interaction between magnetic impurities for two selected graphene nanoflakes containing four hexagonal rings in their structure, differing by their geometry. We describe the electronic structure of either charge-neutral or doped nanoflakes using the tight-binding approximation with the Hubbard term, which is treated within the molecular-field approximation. We find pronounced differences in the RKKY coupling energies, dependent on the placement of the pair of magnetic moments in the nanostructure and on the edge form. For an odd total number of electrons in the structure, we predict in some circumstances the existence of ferromagnetic coupling with leading first-order perturbational contribution, while for an even number of charge carriers the usual, second-order mechanism dominates. Therefore, doping of the nanoflake with a single charge carrier is found to be able to change the coupling from an antiferromagnetic to a ferromagnetic one for some geometries.

###### pacs:

75.30.Hx; 75.75.-c; 73.22.Pr## I Introduction

Two-dimensional graphene, being one of the most promising contemporary materials Geim and Novoselov (2007); Geim (2009) offers unique physical properties due to its electronic structure.Novoselov et al. (2005); Castro Neto et al. (2009); Abergel et al. (2010) Since its discovery it has attracted concerted theoretical and experimental efforts aimed at understanding its rich physics and making basis for future electronics. One of the directions is focused on integrating the charge and spin degrees of freedom, to create novel spintronic devices.Rocha et al. (2010); Yazyev and Katsnelson (2008); Zhang et al. (2011) This encourages studies of magnetic properties of graphene.

Recently, increasing interest has been focused on the graphene nanostructures.Snook and Barnard (2011); Räder et al. (2006); Fernández-Rossier (2008) Such nanostructures can be engineered to merge the properties of graphene with ultrasmall sizes, and can serve as possible building blocks of novel devices (see for example Refs. Ezawa, 2007, 2009a, 2009b; Krompiewski, 2011), especially in context of transport properties, showing giant magnetoresistance.Krompiewski (2009) Significant efforts are made to understand magnetism emerging in finite-size graphene structures without magnetic impurities introduced (an extensive review of that particular topic is presented in the Ref. Yazyev, 2010). However, the possibility of introducing magnetic impurity atoms to the graphene lattice is also taken into account.Santos et al. (2010); Wu et al. (2010); Kang et al. (2011)

Among the topics attracting the interests of researchers, the problem of indirect magnetic coupling between magnetic impurity spins in graphene, mediated by charge carriers, is worth mentioning. This kind of interaction, known as Ruderman-Kittel-Kasuya-Yosida (RKKY) coupling,Ruderman and Kittel (1954); Kasuya (1956); Yosida (1957) can be expected to show unique properties in graphene, different from the behaviour in two-dimensional metals,Béal-Monod (1987) owing to the peculiar, linear dispersion relation for the charge carriers in the vicinity of Dirac points. Various calculations of RKKY interaction in graphene sheets are present in the literature, exploiting the bipartite nature of the graphene crystalline lattice in various (mainly perturbational) approaches.Vozmediano et al. (2005); Dugaev et al. (2006); Cheianov and Fal’ko (2006); Saremi (2007); Bunder and Lin (2009); Brey et al. (2007); Sherafati and Satpathy (2011a, b); Kogan (2011) Moreover, tight-binding calculations in real space were performed for such a system.Black-Schaffer (2010a, b) Recently, also the problem of two Kondo impurities attracted some attention. Uchoa et al. (2011); Hu et al. (2011)

However, the main efforts have so far focused on the calculations of RKKY coupling properties in infinite graphene planes. On the other hand, the coupling in ultrasmall nanoflakes (or nanodisks, Ezawa (2007) containing just a few hexagonal rings) can be expected to deviate from the predictions for infinite system, as the dominance of the edge in a nanoflake substantially modifies the electronic structure and severely breaks the translational symmetry. Especially, the peculiar features of the electronic state at the zigzag edge of graphene (or even graphite) are known both from theoretical predictionsFujita et al. (1996); Nakada et al. (1996) and experimental observations for various systems (e.g. Ref. Klusek et al., 2000), especially quantum-dot systems.Ritter and Lyding (2009) The significance of the edge structure has been already found for example in the transport properties of nanoflakes.Krompiewski (2009) Moreover, in ultrasmall molecule-like structures, the existence of discrete electronic states significantly separated in energy allows us to expect a range of novel phenomena. Therefore, a sound motivation appears for studies of RKKY coupling between magnetic moments in nanosized graphene structures, which is the aim of the present work.

To achieve this goal, we use the real-space tight-binding approximation (TBA) supplemented with the Hubbard term for finite, ultrasmall graphene nanoflakes. We perform the exact diagonalization of the single-particle Hamiltonian resulting from the mean field approximation (MFA). Not limiting the calculations to charge-neutral structures, we take into account the possibility of varying the charge concentration electron by electron (charge doping). One of our goals is to search for the configuration for which the coupling changes its character from ferro- to antiferromagnetic as a result of adding or removing a single electron from the system, in close vicinity to the equilibrium electron concentration. In principle, such a possibility might be opened by placing the nanostructure between two tunneling electrodes, providing a single-electron control of mutual orientation of impurity spins (let us mention that the idea of doping-controlled coupling between magnetic moments in bilayer graphene has been very recently raised in the Ref. Killi et al., 2011). The phenomenon of Coulomb blockade in some graphene nanostructures with two electrodes was studied in Ref. Ezawa, 2008, while another route to modifying the charge concentration is the absorption of gas molecules (as shown in Ref. Schedin et al., 2007).

## Ii Theory

The analysis of RKKY interaction in graphene-based structures consisting of carbon atoms bases on the reliable description of the electronic properties. The most commonly used approach to that problem is the TBA Wallace (1947); Reich et al. (2002); Castro Neto et al. (2009). The main advantage of the TBA method is the ability to include the full, realistic band structure of graphene with finite bandwidth, not being restricted to the linear part of the dispersion relation close to Dirac points. Therefore, the use of certain cut-off schemes is not necessary to obtain convergent results and the issue of the possible unphysical modifications of the RKKY range function is absent. Moreover, the TBA facilitates the non-perturbational treatment of the problem of magnetic impurities. In our study we will adopt this approach, taking into account the electronic hopping for nearest-neighbours only.

The total Hamiltonian for the ultrasmall graphene structure with two magnetic impurities can be written as

(1) |

and contains the contribution from the TBA hopping term , Hubbard term as well as impurity potential .

The TBA term has the following form:

(2) |

Here, and denote creation and annihilation operators for electrons at sites and in the nanoflake, with spin . The summation over nearest-neighbour sites is denoted by . The parameter (usually taken as 2.8 eV) is the hopping integral between nearest neighbours. In the further numerical results, all the energies are normalized to .

In order to incorporate the electron-electron correlations induced by the Coulombic interactions, we include the Hubbard term in the Hamiltonian:

(3) |

with the electron number operators .

The Hubbard term, capturing the on-site coulombic repulsion only, neglects the long-range part of the interaction. The usefulness of this model is a subject of debate (see the recent review in Ref. Kotov et al., 2010). However, it is of noticeably wide use in the theory of carbon (nano)structures; see Refs. Jung and MacDonald, 2009; Yazyev, 2010; Potasz et al., 2011; Fujita et al., 1996; Fernández-Rossier, 2008; Fernández-Rossier and Palacios, 2007 and the recent work on carbon nanotubes in Ref. Alfonsi et al., 2010 Along the lines of the discussion in Ref. Alfonsi et al., 2010, should be treated as an effective parameter describing the influence of Coulombic interactions. Its value might result from a competition between on-site repulsion and the interaction between electrons located on different sites (especially nearest neighbour). Also, some recent results evidence suppression of the influence of the long-range interactions Reed et al. (2010).

Contrary to the choice of hopping integrals in the TBA term, the situation with the value of on-site Coulomb repulsion parameter appears less clear. The values accepted vary between eV Jung and MacDonald (2009) and even eV.Wehling et al. (2011) In our study, if we include the Hubbard term, we accept a moderate value of , close to the choice of Yazyev Yazyev (2010) and Potasz et al. Potasz et al. (2011).

The interaction of the -component of the on-site impurity spin (located in at lattice site ) with an electron spin at the same site is described by the Anderson-Kondo Hamiltonian. For two impurities, located at the sites and , we have:

(4) |

where is a spin-dependent impurity potential (contact potential).

Let us note that we select the Ising form of the interaction Hamiltonian [4] just for simplicity of the calculations. It leads to the final interaction between magnetic impurities described by . The usage of the Heisenberg exchange Hamiltonian, of the form , would just yield , without any modification to the indirect exchange integral itself, which is the only subject of our interest in the present work.

In order to evaluate the RKKY coupling between the impurities at , the ground state energy must be found first in the presence of electrons in the structure, for both parallel and antiparallel orientation of the impurity spins. In undoped graphene each carbon atom donates one electron, so that (i.e. half-filling) characterizes the state of charge neutrality. In general, however, the number of electrons present in the system can vary hypotetically between (empty energetic spectrum) and (completely filled spectrum). We limit our further considerations to , in order not to lose the accuracy of electronic spectrum reproduction by means of the TBA method. The RKKY coupling energy between the impurity spins can be related to the difference of total energies by

(5) |

where the positive value of corresponds to ferromagnetic coupling (F) and the negative value to antiferromagnetic coupling (AF). In our case (when electronic structure description involves NN hopping only) the coupling values calculated for are identical (i.e., electron and hole doping lead to the same results) due to electron-hole transformation symmetry of the Hamiltonian on the finite bipartite lattice.

Let us note that the indirect charge carrier mediated interaction resulting from our calculations can be not necessarily the usual form of RKKY coupling, which is proportional to the square of the contact potential , as resulting from the second order perturbation calculus. Kogan (2011) However, we will use in general the term ”RKKY interaction” to name the indirect charge carrier mediated coupling between impurity spins, even with different characteristic features, resulting from the simultaneous presence of other mechanisms.

In order to deal with the Hubbard term in the Hamiltonian for quite a large system, we adopt the mean field approximation (MFA), which consists in replacement of the form . This approach has been shown recently to compare successfully with some exact diagonalization method for graphene nanostructures Feldner et al. (2010) especially for and has been applied to studies of edge magnetic polarization in graphene sheets Fujita et al. (1996); Jung and MacDonald (2009); Yazyev et al. (2011); Feldner et al. (2011) or RKKY in infinite graphene.Black-Schaffer (2010b) The approximation leads to the effective Hamiltonian defined in single-particle space. The pair of coupled effective single-particle Hamiltonians, and , for spin-up and spin-down electrons, is obtained (the Hamiltonians depend on the electronic densities and can be treated self-consistently).

The Hamiltonian matrix for MFA Hamiltonians can be found in the orthonormal basis of single-electron atomic orbitals . Then, single-particle eigenstates indexed by can be found for each electron spin orientation in the form of linear combination of atomic orbitals . Here, the coefficients for set up an eigenvector of the Hamiltonian matrix corresponding to the eigenvalue .

Let us sort the eigenvalues in an ascending order, so that is the lowest one, etc. Then, the total electronic density at site in the presence of electrons with spin orientation in the system in the ground state at the temperature can be calculated as . The corresponding ground-state energy is

(6) |

The values of and (which add up to the given value of ) should be selected so to minimize the total energy.

After the calculation, the obtained values of electron densities are substituted back into the MFA Hamiltonians and the numerical procedure is repeated iteratively until the satisfactory convergence of the eigenvalues and eigenvectors is achieved. Then the obtained self-consistent numeric value of total energy can be used to calculate RKKY coupling according to the formula given in Eq. 5.

## Iii The numerical results and discussion

For our calculations, we selected two graphene nanoflakes, consisting of four hexagonal rings with an even number of carbon atoms. The first one is similar to a pyrene molecule and contains carbon atoms, and its edge has a zigzag-like character. The second one, triphenylene-like, composed of carbon atoms, possesses the armchair-type edge. Both structures are schematically depicted in Fig. 1(c) and Fig. 7(c).

To study the influence of broken translational symmetry on the RKKY coupling in nanostructures, we calculated between two impurities being either nearest or second neighbours, focusing in each case on two different locations of impurity pairs in a given nanostructure. We denote these cases as 1a,1b and 2a,2b, respectively. The considered pairs of magnetic impurities are presented schematically for both graphene nanoflakes in Fig. 1(c) and Fig. 7(c). Let us mention that nearest neighbours correspond to two ions in two different sublattices while second neighbours constitute a pair of ions in the same sublattice (referring to the subdivision of the bipartite graphene lattice into A and B sublattices for the finite system). In our calculations, we set . In the presentation of further results we focus mainly on the influence of charge doping on the characteristics of indirect coupling.

### iii.1 Pyrene-like nanoflake

In order to gain some insight into the electronic structure of the pyrene-like nanoflake, we plot the predicted energy levels resulting from diagonalization of the Hamiltonian in the absence of magnetic impurities and for in Fig. 1(a). The energy states are numbered by and sorted according to ascending energy and each of the states is doubly (spin) degenerate. No additional degeneracy is observed. In the case of charge neutrality the HOMO-LUMO (highest occupied molecular orbital-lowest unoccupied molecular orbital) gap amounts to 0.89 eV (which is in concert with the value resulting from the calculations in Ref. Feldner et al., 2010, Fig. 2). In order to visualize the electronic densities assigned to the distinct states, we present Fig. 1(b). There, the values of (probabilities of finding the electron at the given lattice site for a given state ) are plotted on the nanoflake scheme, for selected orbitals which are HOMO orbitals for . Let us observe that the subsequent states are characterized by significant variability of the corresponding partial charges distribution. If the selected state is occupied only by a single electron, the distribution of partial charge for this orbital reflects also the spin density.

In Fig. 2 we present the values of RKKY exchange integrals calculated for two positions of magnetic impurities, 1a and 1b, as marked in Fig. 1(c). In order to show qualitatively the influence of increasing contact potential , we prepared the plot for two selected values of and . Moreover, to enable the observation of the Hubbard term effect, we took into consideration the case of as well as . The indirect exchange values are calculated for various deviations in electron number from charge neutrality, for . As can be observed for both impurity pair positions, the coupling is antiferromagnetic for charge-neutral structure (expected for a bipartite lattice half-filled with the electrons, especially for undoped infinite graphene; e.g. Black-Schaffer (2010a)). The interaction is relatively weak and its non-linear rise with increasing values is also observable. The situation is quite similar for , the only difference being the stronger coupling. However, it is quite striking that changing the electronic concentration by one charge carrier from results in switching of the coupling sign from antiferromagnetic to ferromagnetic, which takes place for as well. In these cases, it can be observed that the change of with increasing is noticeably slower than for even values.

In order to study more systematically the dependence of exchange coupling on contact potential , we performed the appropriate calculations for nearest-neighbour ions, the results of which are presented in the Fig. 3 in double logarithmic scale, which allows us to linearize power-law dependences. The case of 1a impurity pair for is illustrated in the Fig. 3(a). It can be verified that two kinds of dependences of on can be distinguished. For selected even total numbers of electrons in the nanoflake () the usual dependence is obeyed, which (according to Fig. 2) coincides with antiferromagnetic sign of interaction. On the contrary, for odd numbers of electrons () we deal with ferromagnetic interaction with . The latter behaviour tends to convert into to the quadratic dependence on (and antiferromagnetic exchange) provided that the potential is strong enough. For clarity, the data for were omitted, as very close to the results for . Let us also observe that qualitatively the same situation is met in presence of the Hubbard term, for [Fig. 3(b)]. It is worth special emphasis that the usual derivation of RKKY exchange integral, within the framework of second-order perturbation calculus, yields always the coupling proportional to the square of the contact potential, which has been recently recapitulated for bipartite graphene lattice. Kogan (2011) Therefore, it is of special interest to identify the source of the unusual linear behaviour. Quite similarly, such a behaviour can be observed for the pair of magnetic impurities in location 1b for example when [Fig. 3(c)] (there, the data for are identical, the same being true for the case of ). Thus, it can be deduced that for selected odd total electron numbers we deal with .

To investigate more deeply the selected case of two impurities in position 1a and and identify the origin of the mentioned unusual behaviour, we plot the energy difference between the AF and F orientations of impurity spins (note that ) as a function of the contact potential in Fig. 4(a). For , the ground state is characterized by unequal number of spin-up and spin-down electrons. In the plot we resolve the contributions coming from the highest energy orbital occupied by a single electron as well as the total contribution originating from the lower energy orbitals occupied by pairs of electrons of opposite spins. It is noticeable that the latter contribution favours the AF state and is proportional to for not too strong contact potentials (as expected on the basis of second-order perturbation calculus). On the contrary, the orbital occupied by a single electron gives rise to the energy difference which is proportional to and tends to prefer the F state of the impurities. In the limit of low , the situation is ruled by the singly coupled state so that the total indirect exchange is ferromagnetic and linear in contact potential, while the increase of leads first to compensation of both contributions and then to domination of the ordinary perturbational -proportional behaviour. The linearity of energy difference in can be explained basing on the fact that for a single electron occupying a given orbital, the leading correction to energy of the orbital [coming from the term Eq. (4) in the Hamiltonian] is of the first- order perturbational kind. If the orbital is occupied by a pair of electrons with opposite spins, first-order corrections for them are of opposite values and cancel each other, while the second-order corrections give rise to the ordinary RKKY interaction. However, when there is no second electron, the uncompensated first-order contribution dominates. Under such circumstances, for F polarization of impurity spins, the first-order correction to the electronic orbital energy due to the presence of two impurities at sites and amounts to (corresponding to the direction of electron spin which minimizes the total energy). Let us assume without loss of generality that . If the polarization of impurity spins is AF, then the corresponding correction is . The resulting energy difference clearly makes a ferromagnetic contribution to the indirect coupling. Its magnitude is proportional to the smaller of the electronic densities on the impurity pair sites; thus it is maximized for equal electronic densities on both sites. Such a contribution to indirect coupling yields some resemblance to the double exchange mechanism. Nolting et al. (2010) However, let us still use the term RKKY interaction to characterize the indirect charge carrier mediated coupling which results from our calculations.

Let us observe, that in some cases we deal with the situation when the electronic density for a given orbital vanishes at least at one of the sites at which the impurity spins are localized. Such an orbital does not indicate any energy difference between AF and F orientation of impurity spins and thus does not give any contribution to the RKKY exchange integral. This is, for example, the case for and impurities in the 1a position; see the corresponding electronic densities for the orbital in Fig. 1(b), which is the highest energy orbital occupied by a single electron. Such a situation prevents the indirect interaction from switching to the F sign (like for ), even though the total number of electrons in the system is odd. Quite a similar situation can be observed for impurities in the 1b position, since the electronic density for the orbital also vanishes at one of the impurity sites. Therefore, the coupling for and is exactly the same. We note that the value of is also unchanged for , which corresponds to the case when the highest energy orbital is occupied by two electrons with opposite spins. However, in that case we can observe that the unperturbed electronic densities are almost equal at both 1b impurity sites. Therefore, the AF state of impurities changes the orbital energy particularly weakly [as can be seen in Fig. 5(b) for ] and the energy changes for spin-up and spin-down electron for F impurity polarization cancel each other. Therefore, a doubly occupied state also gives a particularly weak contribution to exchange integral. This explains the robustness of weak AF RKKY coupling for 1b pair location with respect to deviations from charge neutrality. On the other hand, the same feature of the orbital when it is singly occupied gives rise to a particularly strong ferromagnetic contribution to the coupling (seen clearly for , when ). As mentioned before, this strong ferromagnetic contribution can be explained as a first-order perturbational effect.

To generalize, the possibility of ferromagnetic coupling for nearest-neighbour impurities is open provided that the number of electrons in the system is odd. Under such a condition, the nanoflake gains nonzero magnetic moment, which originates from the HOMO orbital. An additional condition is that the electronic density for the highest energy occupied orbital cannot vanish at both impurity locations. The coupling is particularly enhanced when the electronic densities at the impurity sites are high and close to each other. Then it appears straightforward to explain, that it is energetically favourable to have both impurity spins parallel, as the energy of such a configuration is significantly lowered by the first-order term. Let us notice that in general, for a given weak contact potential , the ferromagnetic couplings, if present, are much stronger than the corresponding antiferromagnetic interactions, which can be seen in Fig. 3.

Let us observe that the presence of a Hubbard term with influences remarkably the interaction, especially for the case of charge neutrality, where it leads to strong enhancement of the AF coupling, and eventually it is able to suppress totally any ferromagnetic behaviour for odd electron number, provided that the contact potential is strong enough. This tendency is observable for both nearest-neighbour impurities in the 1a and 1b positions.

For second-neighbour impurities (results plotted in Fig. 6), it is observed that the interaction for both 2a and 2b impurity locations [see Fig. 1(c)] is almost always ferromagnetic, regardless of the number of electrons, with an exception of . Let us mention that the ferromagnetic coupling is expected for an infinite charge-neutral graphene when considering the impurities belonging to the same sublattice, which is the case for second neighbours (e.g. Ref. Black-Schaffer, 2010a). However, for the nanoflakes, pronounced differences emerge depending on the location of the impurity pair. In the case of the 2a location, the coupling is much enhanced for an odd number of electrons (which can be attributed to the same mechanism as that described for F coupling between nearest neighbours; see the linear dependence of on in Fig. 6 for ). For even values of the interaction energy is much weaker. The presence of the Hubbard term with tends to build up the coupling. The interaction between impurities in the 2b position exhibits traces of vanishing electronic density of the orbital [see Fig. 1(b)], in analogy to the situation met for the 1b pair. Here, the interaction is strongly damped for , while for we observe enhanced F coupling owing to the mechanism mentioned earlier. What is quite interesting, converts the coupling into an antiferromagnetic one.

In order to comment on the influence of the Hubbard term, we can conclude that its dominant result is to enhance the magnitude of the coupling. The Hubbard term associated energy is lowered when the electronic densities for opposite-spin electrons tend to become unequal, hence acting in hand with the impurity potential in creating an imbalance in spin-up and spin-down electronic densities. This could qualitatively justify why the presence of the Hubbard term increases the magnitude of RKKY coupling, as observed for infinite graphene by Black-Schaffer. Black-Schaffer (2010b) Let us exemplify this for the particular case of the 2b impurity pair. The energy differences between the AF and F states of on-site impurities are presented in Figs. 4(b) and 4(c), in the absence and in the presence of Hubbard term, respectively. The total energy difference is plotted together with the separated contribution of the orbitals , each occupied by a pair of opposite-spin electrons. The orbitals are characterized by especially high electronic densities on the impurity sites. Here, all the terms are proportional to , i.e., present an ordinary perturbational behaviour for low . As is visible, the energy differences are much more pronounced for all orbitals in the presence of .

### iii.2 Triphenylene-like nanoflake

The electronic structure of the triphenylene-like nanoflake (obtained from diagonalization of the Hamiltonian in the absence of magnetic impurities and for ) is plotted in Fig. 7(a). What can be observed is that some states acquire additional two-fold degeneracy, in addition to ordinary spin-degeneracy. For charge-neutral nanoflake we have the HOMO-LUMO gap of 1.37 eV. The electronic densities assigned to the distinct states are presented in Fig. 1(b) for selected orbitals which constitute HOMO orbitals when . There, the values of (partial charges present at the lattice sites for a given state) are plotted on the nanoflake scheme. If the selected state is occupied only by a single electron, the distribution of partial charge for this orbital reflects also the spin density.

Figures 8(a) and 8(b) show the values of the RKKY coupling between the impurities in two nearest-neighbour positions 1a and 1b [see Fig. 7(c)]. For an undoped case, an antiferromagnetic interaction is always predicted, again in agreement with expectations for the bipartite lattice. The interesting situation arises for the impurity pair 1b, where the doping with switches the interaction sign to the ferromagnetic one. This sign remains robust against further doping, up to , and also the magnitude of coupling varies only weakly. Moreover, the linear dependence of F interaction on is visible. This phenomenon, present for odd values, sparks off from the mechanism described earlier, involving orbitals occupied with a single electron. Here, it is visible in Fig. 7(b) that the degenerate orbital is characterized by a large electronic density on both 1a impurity sites. Actually, breaking the symmetry of the nanoflake by introducing the magnetic impurities causes that it is most energetically favourable to have large electronic density at both sites of the 1a pair. What is interesting is that the linear dependence of on survives also for . Here, it can be attributed to the fact that for the F polarization of impurity spins, the ground state of the system involves both and orbitals, each occupied by a single electron, instead of just one doubly occupied orbital. Let us note that the electronic densities for degenerate orbitals fulfill the condition , while the orbitals are orthogonal. Therefore, if for one of the orbitals the maximum electronic density at 1a pair sites is achieved, then the second orbital is characterized by vanishing electron density at the same sites. As a consequence, in the F state of the impurities, only one of the orbitals gives contribution to the indirect coupling. The situation is thus analogous to the one expected for the total odd number of electrons. The described occasion of having F coupling for even is intimately connected with an additional degeneracy of some electronic orbitals.

For the 1b pair, the coupling switches its sign from AF to F each time when changes from odd to even. The coupling magnitudes are relatively weak. Ferromagnetic couplings are still linearly dependent on . This time, the degeneracy of the orbital does not play a role at and the coupling is AF.

Let us observe that for the magnetic impurities in the 1a position, the presence of the Hubbard term with tends to enhance the coupling magnitude, conserving its sign. On the other hand, for the 1b pair, the Hubbard term strongly pushes the interaction toward antiferromagnetic for low .

The RKKY exchange integrals calculated for second-neighbour magnetic impurities [in position 2a, 2b; see Fig. 7(c)] are presented in Fig. 9. For 2a placement of the impurities, the interaction is almost always ferromagnetic, but its magnitude is significantly enhanced for odd values of . For the 2b pair of impurities, the coupling is particularly strong (and ferromagnetic) for , which can be attributed to the same mechanism as described previously in the case of the 1b pair.

## Iv Concluding remarks

We have performed a tight-binding based study of RKKY interaction in two graphene nanoflakes, different in their geometries but both containing four hexagonal rings. We focused on the influence of possible charge doping on the magnitude and sign of an indirect interaction. We also incorporated a Hubbard term in the Hamiltonian to estimate the effect of prototypic Coulomb interaction. In general, we found a pronounced dependence of the RKKY coupling integrals on the location of the pair of on-site magnetic moments in the nanostructure.

In relation to the odd number of charge carriers in the nanoflakes, we observed a specific contribution to indirect interaction, originating from the highest-energy electronic orbital occupied by a single electron and giving rise to nonzero spin distributed over the structure. This contribution is proportional to the magnitude of the contact potential and results from an uncompensated first-order correction to the orbital energy due to the presence of a pair of ferromagnetically polarized impurities. This mechanism, which may be regarded as somewhat similar to double exchange, always produces ferromagnetic contribution to indirect exchange and for sufficiently low it can dominate over the typical contribution proportional to originating from second-order perturbational correction to the energy (being the clue of the ordinary RKKY interaction). Such a ferromagnetic contribution can either change the coupling sign from AF to F for nearest-neighbour impurities or enhance the ferromagnetic coupling between second neighbours, as shown for two types of nanoflakes. The details of the mechanism depend also on the presence or absence of degeneracies in energy spectrum of the nanoflake.

The mechanism leading to the indirect ferromagnetic coupling linearly proportional to might be potentially useful, since it allows switching from AF coupling characteristic of nearest-neighbour impurities for charge neutrality to F coupling by adding or removing a single electron from the nanoflake. Especially advantageous conditions for such a situation occur for a triphenylene-like nanoflake for impurities in the 1a position. What is more, adding one or two further charge carriers does not alter the coupling sign in these particular situations.

The influence of the Hubbard term is mainly to enhance the magnitude of the coupling, as noted in Ref. Black-Schaffer, 2010b for infinite graphene. Unequal electronic densities for opposite-spin electrons lower the energy coming from the Hubbard term, and the presence of the local spin-dependent impurity potential acts in the same direction. Therefore, the Hubbard term tends to increase the interaction energy in general. On the other hand, the presence of additional long-range repulsion between the electrons residing on different sites should qualitatively tend to compensate the effect of on-site repulsion. This effect could be captured in treating the parameter not necessarily as a true on-site energy, but rather in the spirit of the effective parameter (as mentioned before, in connection with Ref. Alfonsi et al., 2010)). Sometimes other changes are detected, which can be attributed to the fact that the presence of the Hubbard term itself causes some redistribution of the charge densities over lattice sites with respect to the case when the Coulombic correlations are neglected. These changes also modify the indirect coupling. The application of the MFA approximation to the Hubbard model appears justified for the ground state. Nolting et al. (2010) However, development of the more accurate approaches (albeit consuming less resources than exact diagonalization) would be valuable with a view to the studies of temperature properties of the nanoflakes.

The present calculations are performed for , exploiting ground-state properties of the system. However, in the ultrasmall, molecule-like structures, the separation of the discrete electronic energy levels is usually of the order of tenths of , where eV [see Fig. 1(c) and 7(c)]. Therefore, it appears that no significant redistribution of the charge carriers between the states can happen due to thermal excitations up to the temperatures of interest (i.e., room temperature). As a consequence, the results for coupling energies appear valid also for nonzero temperatures. Perhaps some finite temperature modifications can be expected in the presence of degenerate spectrum, as the impurity-induced energy splitting between the otherwise degenerate states is considerably small (note the interesting thermodynamics of the undoped metallic-like nanoflakes with degenerate zero-energy states, discussed by Ezawa Ezawa, 2007).

In general, we find an immense influence of the electronic structure of the nanostructures on the properties of RKKY interaction, being dependent on the wavefunctions of single orbitals, which are very strongly shaped by the nanoflake geometry. This distinguishes our case from the case of an infinite system, where propagating states are involved, and might open the door for designing the appropriate structures to guarantee the expected indirect coupling features.

###### Acknowledgements.

The author is deeply grateful to T. Balcerzak for critical reading of the manuscript. The numerical calculations have been performed partly with the Wolfram Mathematica 8.0.1 softwareWolfram Research, Inc. (2010) and partly using the LAPACK packageAnderson et al. (1999). This work has been supported by Polish Ministry of Science and Higher Education by a special purpose grant to fund the research and development activities and tasks associated with them, serving the development of young scientists and doctoral students.## References

- Geim and Novoselov (2007) A. K. Geim and K. S. Novoselov, Nat. Mater. 6, 183 (2007).
- Geim (2009) A. K. Geim, Science 324, 1530 (2009).
- Novoselov et al. (2005) K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, M. I. Katsnelson, I. V. Grigorieva, S. V. Dubonos, and A. A. Firsov, Nature (London) 438, 197 (2005).
- Castro Neto et al. (2009) A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, Rev. Mod. Phys. 81, 109 (2009).
- Abergel et al. (2010) D. S. L. Abergel, V. Apalkov, J. Berashevich, K. Ziegler, and T. Chakraborty, Adv. Phys. 59 (2010).
- Rocha et al. (2010) A. R. Rocha, T. B. Martins, A. Fazzio, and A. J. R. da Silva, Nanotechnology 21, 345202 (2010).
- Yazyev and Katsnelson (2008) O. V. Yazyev and M. I. Katsnelson, Phys. Rev. Lett. 100, 047209 (2008).
- Zhang et al. (2011) Q. Zhang, K. S. Chan, and Z. Lin, Appl. Phys. Lett. 98, 032106 (2011).
- Snook and Barnard (2011) I. Snook and A. Barnard, Physics and Applications of Graphene - Theory, edited by S. Mikhailov (InTech, Rijeka, 2011), Chap. 13, p. 277.
- Räder et al. (2006) H. J. Räder, A. Rouhanipou, A. M. Talarico, V. Palermo, P. Samori, and K. Müllen, Nat. Mater. 5, 276 (2006).
- Fernández-Rossier (2008) J. Fernández-Rossier, Phys. Rev. B 77, 075430 (2008).
- Ezawa (2007) M. Ezawa, Phys. Rev. B 76, 245415 (2007).
- Ezawa (2009a) M. Ezawa, Eur. J. Phys. B 67, 543 (2009a).
- Ezawa (2009b) M. Ezawa, New J. Phys. 11, 095005 (2009b).
- Krompiewski (2011) S. Krompiewski, Central Eur. J. Phys. 9, 369 (2011).
- Krompiewski (2009) S. Krompiewski, Phys. Rev. B 80, 075433 (2009).
- Yazyev (2010) O. V. Yazyev, Rep. Progr. Phys. 73, 056501 (2010).
- Santos et al. (2010) E. J. G. Santos, D. Sánchez-Portal, and A. Ayuela, Phys. Rev. B 81, 125433 (2010).
- Wu et al. (2010) M. Wu, C. Cao, and J. Z. Jiang, New J. Phys. 12, 063020 (2010).
- Kang et al. (2011) J. Kang, H.-X. Deng, S.-S. Li, and J. Li, J. Phys. Condens. Matter 23, 346001 (2011).
- Ruderman and Kittel (1954) M. A. Ruderman and C. Kittel, Phys. Rev. 96, 99 (1954).
- Kasuya (1956) T. Kasuya, Progr. Theor. Phys. 16, 45 (1956).
- Yosida (1957) K. Yosida, Phys. Rev. 106, 893 (1957).
- Béal-Monod (1987) M. T. Béal-Monod, Phys. Rev. B 36, 8835 (1987).
- Vozmediano et al. (2005) M. A. H. Vozmediano, M. P. López-Sancho, T. Stauber, and F. Guinea, Phys. Rev. B 72, 155121 (2005).
- Dugaev et al. (2006) V. K. Dugaev, V. I. Litvinov, and J. Barnas, Phys. Rev. B 74, 224438 (2006).
- Cheianov and Fal’ko (2006) V. V. Cheianov and V. I. Fal’ko, Phys. Rev. Lett. 97, 226801 (2006).
- Saremi (2007) S. Saremi, Phys. Rev. B 76, 184430 (2007).
- Bunder and Lin (2009) J. E. Bunder and H.-H. Lin, Phys. Rev. B 80, 153414 (2009).
- Brey et al. (2007) L. Brey, H. A. Fertig, and S. Das Sarma, Phys. Rev. Lett. 99, 116802 (2007).
- Sherafati and Satpathy (2011a) M. Sherafati and S. Satpathy, Phys. Rev. B 83, 165425 (2011a).
- Sherafati and Satpathy (2011b) M. Sherafati and S. Satpathy, Phys. Rev. B 84, 125416 (2011b).
- Kogan (2011) E. Kogan, Phys. Rev. B 84, 115119 (2011).
- Black-Schaffer (2010a) A. M. Black-Schaffer, Phys. Rev. B 81, 205416 (2010a).
- Black-Schaffer (2010b) A. M. Black-Schaffer, Phys. Rev. B 82, 073409 (2010b).
- Uchoa et al. (2011) B. Uchoa, T. G. Rappoport, and A. H. Castro Neto, Phys. Rev. Lett. 106, 016801 (2011).
- Hu et al. (2011) F. M. Hu, T. Ma, H.-Q. Lin, and J.-E. Gubernatis, Phys. Rev. B 84, 075414 (2011).
- Fujita et al. (1996) M. Fujita, K. Wakabayashi, K. Nakada, and K. Kusakabe, J. Phys. Soc. Jpn. 65, 1920 (1996).
- Nakada et al. (1996) K. Nakada, M. Fujita, G. Dresselhaus, and M. S. Dresselhaus, Phys. Rev. B 54, 17954 (1996).
- Klusek et al. (2000) Z. Klusek, Z. Waqar, E. Denisov, T. Kompaniets, I. Makarenko, A. Titkov, and A. Bhatti, Appl. Surf. Sci. 161, 508 (2000).
- Ritter and Lyding (2009) K. A. Ritter and J. W. Lyding, Nat. Mater. 8, 235 (2009).
- Killi et al. (2011) M. Killi, D. Heidarian, and A. Paramekanti, New J. Phys. 13, 053043 (2011).
- Ezawa (2008) M. Ezawa, Phys. Rev. B 77, 155411 (2008).
- Schedin et al. (2007) F. Schedin, A. K. Geim, S. V. Morozov, E. W. Hill, P. Blake, M. I. Katsnelson, and K. S. Novoselov, Nat. Mater. 6, 652 (2007).
- Wallace (1947) P. R. Wallace, Phys. Rev. 71, 622 (1947).
- Reich et al. (2002) S. Reich, J. Maultzsch, C. Thomsen, and P. Ordejón, Phys. Rev. B 66, 035412 (2002).
- Kotov et al. (2010) V. N. Kotov, B. Uchoa, V. M. Pereira, A. H. Castro Neto, and F. Guinea (2010), arXiv:1012.3484v1, to appear in Rev. Mod. Phys.
- Jung and MacDonald (2009) J. Jung and A. H. MacDonald, Phys. Rev. B 79, 235433 (2009).
- Potasz et al. (2011) P. Potasz, A. D. Güçlü, O. Voznyy, J. A. Folk, and P. Hawrylak, Phys. Rev. B 83, 174441 (2011).
- Fernández-Rossier and Palacios (2007) J. Fernández-Rossier and J. J. Palacios, Phys. Rev. Lett. 99, 177204 (2007).
- Alfonsi et al. (2010) J. Alfonsi, G. Lanzani, and M. Meneghetti, New J. Phys. 12, 083009 (2010).
- Reed et al. (2010) J. P. Reed, B. Uchoa, Y. I. Joe, Y. Gan, D. Casa, E. Fradkin, and P. Abbamonte, Science 330, 805 (2010).
- Wehling et al. (2011) T. O. Wehling, E. Şaşıoğlu, C. Friedrich, A. I. Lichtenstein, M. I. Katsnelson, and S. Blügel, Phys. Rev. Lett. 106, 236805 (2011).
- Feldner et al. (2010) H. Feldner, Z. Y. Meng, A. Honecker, D. Cabra, S. Wessel, and F. F. Assaad, Phys. Rev. B 81, 115416 (2010).
- Yazyev et al. (2011) O. V. Yazyev, R. B. Capaz, and S. G. Louie, Phys. Rev. B 84, 115406 (2011).
- Feldner et al. (2011) H. Feldner, Z. Y. Meng, T. C. Lang, F. F. Assaad, S. Wessel, and A. Honecker, Phys. Rev. Lett. 106, 226401 (2011).
- Nolting et al. (2010) W. Nolting, and A. Ramakanth, Quantum Theory of Magnetism (Springer, Berlin, Heidelberg, 2009).
- Ezawa (2007) M. Ezawa, Phys. Rev. B 79, 241407(R) (2009).
- Wolfram Research, Inc. (2010) Wolfram Research, Inc., Mathematica, Version 8.0 (Wolfram Research, Inc., Champaign, IL, 2010).
- Anderson et al. (1999) E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney, et al., LAPACK Users’ Guide (Society for Industrial and Applied Mathematics, Philadelphia, PA, 1999), 3rd ed., ISBN 0-89871-447-8 (paperback).