A first-principles study on magnetic coupling between carbon adatoms on graphene

The weak ferromagnetism reported for graphite and related carbon nanostructures is frequently related to the magnetic coupling of point defects such as vacancies or hydrogen adatoms that interact with only one sublattice in the bipartite graphene lattice. In this paper, using density functional theory calculations we study the magnetic coupling between point defects, such as carbon adatoms, which form bonds with atoms located on both sublattices. We show that there is ferromagnetic coupling for small separations between the adatoms. Further, we demonstrate that it is energetically favorable for C adatoms to agglomerate. Our results indicate that the magnetism induced in graphite by irradiation can be explained as being due to the existence of adatom agglomerations and small graphene flakes formed from isolated adatoms rather than uniformly distributed point defects.


Introduction
Numerous theoretical predictions (for an overview, see [1,2]) of intrinsic magnetism in carbon materials such as graphite, fullerenes, nanotubes, nanodiamonds and graphene have stimulated a considerable body of experimental work on this subject. Indeed, the experimental confirmation of magnetism in low-density biocompatible carbon systems not involving any d-elements (but possibly for hydrogen and nitrogen impurities), which have traditionally been believed to show diamagnetic behavior, would be very important not only in the context of the fundamental solid state physics, but also from the viewpoint of practical applications.
However, although there have been several reports [3]- [10] on magnetism in carbon materials, the results are quite controversial and not always reproducible, as has been already pointed out [1]. There is very little solid experimental evidence for the existence of ferromagnetism in carbon systems (such as polymerized fullerenes [11,12], irradiated graphite by protons [13,14] or heavy ions [15]- [17]). Many of the results appear to be controversial, as they could not be reproduced by other groups; see [1] for details.
Most previously published papers argue that magnetism originates from defects [18,19] in the graphitic network, which give rise to local magnetic moments. Examples are vacancies [20]- [26], interstitials [27], carbon adatoms [17,28] and atoms on the edges of graphitic nano-fragments with dangling bonds either passivated with hydrogen atoms [29]- [31] or free [30,32]. The nuclear magnetic resonance experiments indeed provide evidence [18] that defects in carbon systems have magnetic moments, as they couple to implanted Fe atoms. Thus, a key issue is the mechanism of the coupling of magnetic moments localized at defects.
The problem of defect coupling has been addressed in several papers [23]- [25] in the framework of density functional theory (DFT). All of them dealt with point defects that affect a single site in the graphenic network, such as vacancies or H adatoms. It was shown that defects that reside on the same sublattice may couple ferromagnetically, while antiferromagnetic (AFM) coupling is more favorable for defects residing on different sublattices [33,34], and that coupling may occur through the Ruderman-Kittel-Kasuya-Yoshida (RKKY) interaction. However, one has to keep in mind that defects may be randomly distributed over the two sublattices, which should give rise to zero net magnetic moment, according to the Lieb theorem [35]. In addition, the RKKY interaction is expected to be weak due to the semi-metallic electronic structure of graphene [33], [36]- [38]. 3 At the same time, the macroscopic ferromagnetic (FM) state may come from long-range interaction between other defects that affect both sublattices, such as carbon and possibly other adatoms on graphene that occupy the bridge [28,39] or hollow (middle-of-the-hexagon) [40] positions. However, the coupling of such defects has never been studied by first-principles methods, but only by the tight-binding method with a model Hamiltonian, which ignores relaxations of the atomic structures [41].
In this paper, in order to assess the magnetic coupling between defects that affect both sublattices in the bipartite lattice of graphene, we carried out extensive computer atomistic simulations based on the DFT approach. We concentrate on carbon adatoms on a graphene sheet as an example of such defects. By evaluating the difference between the FM and the corresponding AFM state, we show that there is FM coupling for small separations between the adatoms. We further demonstrate that it is energetically favorable for C adatoms to agglomerate. However, in some specific configurations, strong rehybridizations of sp 2 and p z orbitals destroy FM ordering.
This paper is organized as follows. Details of our first-principles calculations are given in section 2, and the results of calculations are presented in section 3, where we analyze the adsorption of two carbon adatoms in various configurations and the existence of a coupling between localized magnetic moments. The summery and conclusion are given in section 4.

Computational details and models
We have carried out DFT calculations using the Vienna ab initio simulation package VASP [42,43]. The code uses the full-potential projector augmented waves (PAW) framework [44,45]. The Perdew-Burke-Ernzerhof (PBE) [46] functional, a generalized gradient approximation (GGA)-type functional, was used in order to account for exchange-correlation effects in our spin-polarized calculations. A kinetic energy cutoff of 400 eV was found to be sufficient for achieving a total energy convergence within several meV considering k-point sampling. In the (7 × 7) simulation cell, for (7 × 7) primitive graphene cells, used in most of our calculations, a (3 × 3 × 1) k-point grid was used. We checked our results for several configurations using larger meshes, up to (7 × 7 × 1), and did not find considerable differences in the coupling energies.
All atoms were fully relaxed until the forces on individual atoms were smaller than 0.02 eV A −1 . Considering the smearing issue, during the geometry optimization runs a Gaussian smearing with a width of 0.2 eV was used, while the tetrahedron method with Blöchl corrections [47] was used for accurate energy calculations. The diffusion barriers were estimated by the climbing image nudge elastic band method (CI NEB) [48,49], with a force tolerance of 0.02 eV A −1 and five intermediate geometries for the transition state search.
To find a possible proof of magnetic in-plane interactions between defects in carboneous nano-structures, we have focused our attention on a single graphene layer. Our model can be directly used for explaining experimental results obtained in graphene and with a certain caution to multi-layer systems, as recent theoretical works indicate the importance of inter-layer interaction to magnetic effects [50,51].
We have used a slab supercell calculation with 18.3 Å of vacuum in the c direction as a model of the graphene layer. It ensures absence of reciprocal interaction between periodic images in this particular direction. Note that pristine graphene consists of two equivalent sublattices, denoted A and B in the following, whereas in graphite the sublattices are nonequivalent due to inter-layer interactions [24].  Table 1. Structural, energetic and magnetic properties in As configuration. d ad is the distance between the two adatoms, d 1st is the distance between the µ 2 atom and its first neighbor, d b is the distance between the two basal C atoms bridged by the adatom, µ tot is the total magnetic moment of the cell and µ ad is the projected magnetic moment of the adatoms. Taking into account that the interaction between adatoms could be significantly influenced by the effects of periodic images of the cell, a set of test calculations has been performed by varying the cell size and the k-point sampling accordingly. For instance, in the standard (7 × 7) cell used for the section 3 discussion, a (3 × 3 × 1) k-point grid was used. We have restricted our test case to a specific configuration of two adsorbed C atoms separated by two hexagon-widths, denoted 'As configuration' in the following and depicted in figure 1. Table 1 summarizes the main geometric, energetic and magnetic properties of the particular As configuration in the FM state for several cell sizes. The smallest cell (4 × 4) contains 34 atoms (32 atoms in the graphene sheet and two adatoms), while the largest (10 × 10) cell has 202 atoms. Note that the distances between the atoms are given in the FM state only. However, very small deviations (less than 0.02 Å ) have been observed in the AFM state geometries. As shown previously [28], [52]- [54], a C adatom bridges two carbon atoms of different sublattices, resulting in a µ 2 position, with the perpendicular distance from the adatom to the graphenic surface of 1.87 Å. Our calculations for a single C adatom performed with the PBE functional give almost the same distance, i.e. 1.84 Å. Moreover, the distance between the adsorbed carbon atom and its first neighbor (d 1st ) is about 1.52 Å, whereas the distance between the two basal C atoms (d b ) belonging to A and B sublattices is 1.57 Å. The total magnetic moment of the cell in the case of a single carbon adatom is 0.42µ B , whereas the value in the literature is about 0.45µ B [28,54]. In figure 2, the effects of the adsorption of a single C atom are depicted in the spin-polarized band structures. In addition to the band splitting due to symmetry breaking, two interesting features can be extracted from the spin-polarized band structures. Due to the presence of the localized state close to the Fermi level, the graphene layer becomes metallic, with a slight left-shift of the band crossing initially in the K point, lying between the M and K points of the Brillouin zone. Note that the corresponding anti-bonding state in the spin-down polarization spectrum is also well localized and is further away from the Fermi level. The highest occupied state and the lowest unoccupied band are almost flat in the K point, whereas these states are known to have linear dispersion around this special point [55]. Thus, it is likely that C adatoms affect electronic transport in graphene. When another C atom comes close to the previous adsorbed atom, corresponding to the As configuration shown in figure 1, one observes almost no change in d 1st and d b , with a separation of 4.97 Å between the two adatoms, denoted d ad in the following. Interestingly, the influence of cell size on the main geometric parameters, namely d ad , d 1st and d b , is rather small. For instance, the value of d ad is increased by only 0.07 Å when passing from the smallest calculation cell to the largest one. One possible explanation for this difference is the better localization of the adatom's state in the largest cell. Indeed, at the same time, one observes an increase of the total magnetic moment of the cell from 1.03µ B to 1.77µ B , but more importantly in the projected magnetic moment of the adatoms (0.22-0.29µ B ) for the (4 × 4) and (10 × 10) systems, respectively. As a consequence of this better localization, the repulsion between the adatoms is enhanced.
However, from figure 3 it can be seen that the convergence of the total magnetic moment of the cell with respect to the cell size is non-monotonic. Hence the moment seems to be much more dependent on the system size than on any other characteristic i.e. geometry and even energetics. Indeed, when varying the concentration of adatoms while keeping the system in Asposition, the FM state is still the ground state and the energy difference between the FM and AFM states lies around −35 meV since the cell size is larger than (6 × 6) primitive cells. In other words, there is almost no influence of the periodically repeated images of the adatoms if the distance to the next nearest neighboring adatom is larger than 12 Å. Hence the long-range  nature of the interaction between impurities, usually evoked [56,57], becomes questionable as soon as the adsorbed species are shared by the two sublattices. Taking into account these results, in what follows we restrict our attention to the (7 × 7) primitive cell model, which gives a dilution of the defects of about 2%. Note that this value is still at least 10 times as large as recent experimental estimations [4] of defect concentrations.

Two carbon adatoms
To study the influence of the adsorption position of the second adatom in the neighborhood of the first C adatom on the structural, energetic and magnetic characteristics of the system, we have carried out a series of calculations for the configurations depicted in figure 1.

Structural characteristics
From the geometries of the adsorbed atoms described by d 1st and d b , one can immediately show two distinct types of behavior. Speaking in terms of bond lengths, when the two adatoms are separated by a distance larger than 5 Å, that is, As, E, F, G, H and I configurations, the atomic geometry of the system is very close to those obtained in the single adatom case: d 1st is about 1.53 Å, d b is very close to 1.57 Å, while the distance to the graphene plane is slightly larger than 1.90 Å on average. In these configurations, the effects of adsorption on the geometry are erased in a very short range, i.e. after the second neighbor ring.
In contrast, when the separation between the two adatoms is smaller than 5 Å, for A, B, C and D configurations, strong geometry distortions are observed; see figure 4 for A and C configurations. Due to the lateral interactions of the adatoms, they tend to maximize their separation, inducing a higher elevation of the neighboring atoms away from the graphene plane. This is exacerbated in the A configuration, in which the orbitals clearly change from the sp 2 to a more pronounced sp 3 character, resulting in a buckled structure. Note that, in all tested configurations, very small differences in the geometric parameters (less than 0.1 Å on average) of the FM and AFM states have been observed. From a structural point of view, the adatoms do not feel each other (no changes in the bond lengths between the adatom and neighboring atoms occur) as soon as the distance between them is larger than 5 Å. in meV, with respect to the distance between the two adatoms. Circled points stand for AF states energetically favored.

Magnetic coupling versus distances
In all tested configurations, the magnetic states, i.e. FM or AFM states, have lower energy than the non-magnetic (NM) states by at least 20 meV in the case of C configuration, up to 160 meV in the case of A configuration, and thus are more stable. On average, the NM state lies at about 100 meV above the most stable magnetic state. The interaction of the localized moments is evaluated by calculating the energy difference E FM−AFM between the FM state and AFM state energies. In figure 5, E FM−AFM as a function of the distance between the adatoms is shown. Except for the C configuration, which appears to be very specific due to strong rehybridization (it will be discussed in the next subsection), when adatom separation is less than 6 Å, the FM states are the ground states. On average, the energy difference for these configurations is about 50 meV. By relating the energy difference to the Heisenberg model with a coupling value of J = E FM−AFM /2 per pair of spins, one can estimate the critical temperature T c using the mean-field expression T c = JS(S + 1)/3k B . With S = 1, one obtains a T c value of about 200 K. Now, if one considers the formula, T c = 0.44JS(S + 1)/k B , used in [25] that corresponds to a quantum Monte-Carlo study of a quasi-two-dimensional Heisenberg ferromagnet on a square lattice [58], one gets T c = 255 K.
This rough estimate has only the same order of magnitude as the published experimental values lying in the range of 300-500 K; see [3,4,6]. This discrepancy may be attributed to the well-known underestimation of GGA-flavor DFT calculations of the exchange coupling constants. Indeed, in the case of hybrid functional-type calculations for H-decorated vacancies in interaction, at a distance of about 5 Å, J goes up to 150 meV [25]. However, our estimated value of J, i.e. 25 meV, compares relatively well with the values obtained by GGA-type DFT. Although the nature of defects can vary, a rough estimate of the average value is 40 ± 20 meV. For instance, in the work of Xia et al [17], J is 38 meV, whereas in [59], 28 meV was obtained for two metallic adatoms separated by 5 Å, in the top position on a carbon nanotube's wall. Compared with the value obtained in an analytical model, there is fairly good agreement; see for instance [37]. When the distance between impurities is above 6 Å, the AFM state has the same or even lower energy than the FM state, except for the H configuration. However, one should keep in mind that those energy differences (about 5 meV) are close to the precision limit of our calculations. In addition, the second-nearest neighbors (images of the adatoms due to the use of periodic boundary conditions) may also affect the results. Indeed, for the F, G and H configurations in (7 × 7) cell calculations, the distance to the first second-nearest neighbor is 11.2, 9.8 and 10.7 Å, respectively. Nevertheless, for the (10 × 10) primitive cell slab, the distance between second-nearest neighbors increases and is 18.6, 17.3 and 17.8 Å for the F, G and H configurations, respectively, but only minor changes in E FM−AFM are observed, less than a few meV for each tested case. This may be considered as an indirect proof of the destruction of FM ordering for separations above 6 Å. Thus, it is reasonable to state that above 6 Å no magnetic ordering exists, and that ferromagnetism can appear at smaller separations of the adatoms.

C configuration: the antiferromagnetic state stable at short distance?
The previous statement may be referred to as one resulting in a discrepancy: the C configuration. Indeed, at this specific geometry the minimization of the system energy always gives the AFM state, independent of whether the starting configuration corresponds to the FM or AFM state. Adatom magnetic moments are ±0.16µ B . In figure 4, the side view of the optimized C geometry is given in comparison to the A relaxed structure. Due to the presence of the adatoms that repel each other in the A configuration, the elevation, with respect to the graphene plane, of two neighboring atoms bound to each adatom is relatively high, i.e. 0.9 Å. As a consequence, these specific atoms possess a more pronounced sp 3 character, which explains the drastic change in band structure (see figure 7) when compared to the single adatom case. One important feature of the band structure in the A geometry is the metallic character of the system, with the highest occupied state crossing the Fermi level between and M points and K and points. The second main characteristic is the presence of a flat band very close to the Fermi level in the spin-up density, corresponding to a localized state on the adatoms, responsible for the total magnetic moment of the system. Note that due to this particular state, band crossing in the K points is no longer observed. In summary, for the A configuration the change in hybridization of basal atoms (bound to the adatoms) stabilizes the spin-up states and ensures the FM state. When comparing the geometry of the C configuration with that of A, the main difference is the distance of the basal atoms from the graphene layer, which is decreased by 0.3 Å. As a result, the basal atoms have a less pronounced sp 3 character, which destabilizes the localized states. The extra electrons are shared over several atoms, and the associated bands end up with a gap (0.27 eV) at the M point; see the last panel of figure 8. The highest occupied state is more delocalized, and its dispersion is more pronounced than in the B configuration. Indeed, in B geometry, occupied spin-up states with the two highest energies are very flat and the highest occupied state outcrops the Fermi level at the K point. That also explains the large value of the total magnetic moment obtained in the B configuration (2.0 versus 1.70µ B for A configuration).
However, the C configuration is less stable than the A configuration, with an energy difference of 225 meV. Moreover, when one estimates the energy barrier associated with the adatom migration process, with one adatom kept fixed, from the C to A configuration, one obtains only 30 meV. Comparatively, when passing from B to A configuration, again with one adatom in a fixed position, the energy barrier is larger: about 250 meV. As a consequence, the A configuration is thermodynamically a more probable structure than C, and from the kinetic point of view, the existence of the AFM state in the C configuration is highly improbable due to the high mobility of the adatoms.

Agglomerations of the adatoms
After having estimated the absolute energy of the most stable state, be it the AFM or the FM state, for all tested configurations, it appears that A-position has the lowest energy; see In contrast, positions As, F, G, H and I are less stable than a state with two well-separated adatoms. In the most favorable case, the stabilization energy is about 125 meV per adatom, whereas for the less stable state it is close to 75 meV per adatom. Interestingly enough, all configurations having distances between adatoms of less than 5 Å are more stable than infinitely separated adatoms. This may be considered as evidence for the existence of a driving force to agglomerate adatoms in a small region in the graphitic network, with the impurities shared by the two sublattices. However, a small increase in energy at separations of 6-9 Å may prevent adatom clustering. To study this issue in more detail, we calculated the energy barrier of adatom diffusion, which proved to be 0.52 eV, in good agreement with [28] for the hopping process of a single C adatom. We found that when a second adatom is present in the vicinity of the moving adatom, for instance when passing from B to A configuration, the energy barrier associated with the migration process is almost divided by a factor of two, i.e. 0.25 eV. Thus, provided that the concentration of adatoms is relatively high (otherwise the probability for adatoms to 'find' each other is too low), clustering of adatoms is possible. Keeping in mind the ion impact as a source of adatoms, another necessary condition is the proper temperature range (close to room temperature): if the temperature is too low, adatoms do not move, and if it is too high, adatoms will quickly diffuse away from the local area where they have been created by the energetic ion.
Ion irradiation experiments carried out at low temperatures followed by in situ magnetic measurements at various temperatures may be used to create magnetic defects in graphene and detect them. Then, by increasing the temperature and measuring the magnetic response of the system, the contribution of different types of defects and their agglomeration can be probed. For graphene flakes deposited on TEM grids, the defects can be created and even directly observed [60]. Adatoms should contribute to magnetism at low temperatures when they are immobile. In practice, magnetic effects may come from the co-existence of several types of defect, e.g. adatoms and vacancies. Advanced techniques such as spin-polarized scanning-tunneling spectroscopy may be used to detect the magnetic moments of adatom agglomerations and shed light on their interactions. In particular, superparamagnetism may be observed if the adatom agglomerations are weakly coupled. Moreover, different Curie temperatures may be found provided that agglomerations of adatoms are large enough.

Conclusion
We carried out extensive first-principles computer atomistic simulations aimed at evaluating the coupling between point defects on the graphene lattice and their energetics. In contrast to previous studies that considered defects such as vacancies or hydrogen adatoms that interact with only one sublattice in the bipartite graphene lattice, we studied the magnetic coupling between point defects that form bonds with atoms located on both sublattices in graphene. We concentrated on carbon adatoms on a graphene sheet as an example of such defects. By evaluating the difference between the FM and the corresponding AFM state, we showed that there is FM coupling for small separations between the adatoms. Although C adatoms on graphene sheets may be too mobile to exist at room temperature, such configurations may be stabilized in graphite or graphene paper, a recently manufactured analogue of nanotube buckypaper, which shows interesting magnetic properties [61]. Further, we demonstrated that it is energetically favorable for C adatoms to agglomerate over small distances. On the basis of this observation, one can assume that if the concentration of adatoms is high or the temperature is relatively low, a magnetic island exhibiting intriguing magnetic behavior, e.g. superparamagnetism, can be formed. Adatoms may coalesce, so small graphene flakes between graphite layers (or on top of a graphene sheet) may eventually be formed. Such flakes also exhibit magnetic properties [62], with the spin value being dependent on their shapes due to topological frustration. Overall, our results indicate that the magnetism observed in graphite can be explained as being due to the existence of local defect agglomerations and small graphene flakes with finite magnetic moments rather than being due to uniformly distributed point defects.
work was also performed using HPC resources from GENCI-[CINES] (grant no. 2010-6005). This work was supported by the Academy of Finland through several projects.