Inter-electron interactions and the RKKY potential between H adatoms in graphene

TL;DR

This study uses advanced simulations to explore how electron interactions influence the magnetic coupling between hydrogen adatoms on graphene, revealing effects on potential barriers, ordering, and lattice stability.

Contribution

It provides new insights into how inter-electron interactions modify RKKY interactions and the stability of adatom arrangements on graphene.

Findings

Inter-electron interactions significantly alter RKKY potentials.

Anti-ferromagnetic and charge density wave orders impact RKKY differently.

Certain adatom superlattices show stability against small displacements.

Abstract

We use first-principles Quantum Monte-Carlo simulations to study the Ruderman-Kittel-Kasuya-Yosida (RKKY) interaction between hydrogen adatoms attached to a graphene sheet. We find that the pairwise RKKY interactions at distances of a few lattice spacings are strongly affected by inter-electron interactions, in particular, the potential barrier between widely separated adatoms and the dimer configuration becomes wider and thus harder to penetrate. We also point out that anti-ferrromagnetic and charge density wave orderings have very different effects on the RKKY interaction. Finally, we analyze the stability of several regular adatom superlattices with respect to small displacements of a single adatom, distinguishing the cases of adatoms which populate either both or only one sublattice of the graphene lattice.

Equations

$$\hat{H}=\sum _{⟨x,y⟩,\sigma }-t_{xy}\left({\hat{a}}_{y,\sigma }^{†}{\hat{a}}_{x,\sigma }+h.c.\right)+\frac{1}{2}\sum _{x,y}{V}_{xy}{\hat{q}}_x{\hat{q}}_y,$$

$$(x_1+L_1,x_2)\to (x_1,x_2),$$

$$(x_1,x_2+L_2)\to (x_1+L_2/2,x_2),$$

$${\hat{H}}_{hyb.}=\gamma \sum _{x\in \mathcal{H},\sigma }\left({\hat{a}}_{x,\sigma }^{†}{\hat{c}}_{x,\sigma }+h.c.\right)+E_d\sum _{x\in \mathcal{H},\sigma }{\hat{c}}_{x,\sigma }^{†}{\hat{c}}_{x,\sigma },$$

$$\mathcal{Z}=\int \mathcal{D}{\varphi }_{x,\tau }{e}^{-S\left[{\varphi }_{x,\tau }\right]}\det \left(M_e\left[{\varphi }_{x,\tau }\right]\right)\det \left(M_h\left[{\varphi }_{x,\tau }\right]\right),$$

$$\det \left(M_e\left[{\varphi }_{x,\tau }\right]\right)\det \left(M_h\left[{\varphi }_{x,\tau }\right]\right)={|\det \left(M\left[{\varphi }_{x,\tau }\right]\right)|}^2$$

$${\mathcal{F}}_{xy}=-T\sum _{{ϵ}_{xy}}\ln \left(1+e^{-{ϵ}_{xy}/T}\right)+F_0.$$

$$\Delta \mathcal{F}=-T{\int }_0^{1}d\alpha {\partial }_{\alpha }\log {\mathcal{Z}}_{\alpha },$$

$${\mathcal{Z}}_{\alpha }=\int \mathcal{D}{\varphi }_{x,\tau }{e}^{-S\left[{\varphi }_{x,\tau }\right]}|\det \left(M_{\alpha }\left[{\varphi }_{x,\tau }\right]\right){|}^2,$$

$$\Delta \mathcal{F}=-2T{\int }_0^{1}d\alpha ⟨\mathrm{Re}\mathrm{Tr}\left(M_{\alpha }^{-1}{\partial }_{\alpha }{M}_{\alpha }\right)⟩,$$

Full text

I Introduction

Functionalization of graphene with hydrogen adatoms or other admolecules which produce resonant scattering centers is currently a subject of intense research. First of all, it provides a way to create a band gap in graphene Abanin et al. (2010); Putz et al. (2014); Ulybyshev and Katsnelson (2015) with a possibility to tune it and even to return the material to the initial semi-metallic state Elias et al. (2009). Also the magnetic moments induced around hydrogen adatoms due to inter-electron interactions Yazyev and Helm (2007); Wang and Pantelides (2012); Ulybyshev and Katsnelson (2015) play an important role in spin relaxation processes Han et al. (2014) and can be used to tune the magnetic properties of graphene.

The spatial distribution of adatoms plays a crucial role in the properties of hydrogenated graphene. For instance, magnetic moments of adatoms placed at different sublattices are coupled anti-ferromagnetically, while adatoms placed sufficiently close to each other on the same sublattice induce ferromagnetic ordering Ulybyshev and Katsnelson (2015). The stability (or instability) of these adatom configurations determines the magnetic properties of the functionalized material and might also explain the arrangement of hydrogen adatoms in regular superlattices observed in recent experiments Balog et al. (2010); Lin et al. (2015). Especially important is the case of functionalized graphene on top of boron nitride Balog et al. (2010), where hydrogen adatoms tend to occupy only one sublattice at special places in a moiré structure, forming islands of graphane.

The smallness of the pairwise elastic interaction of hydrogen adatoms in graphene, which does not exceed $\sim 10\,{\rm meV}$ for distances larger than the inter-atomic spacing Branicio et al. (2016) and decays as $r^{-3}$ at large distances Peyla and Misbah (2003), suggests that the Ruderman-Kittel-Kasuya-Yosida (RKKY) contribution from conduction electrons might dominate the inter-adatom interactions in graphene. The RKKY potential between pairs of adatoms was studied in a number of papers starting from the seminal article Saremi (2007). Typically, some analytic approximations to the fermionic Green’s function in presence of resonant scatterers are used Shytov et al. (2009); Lawlor et al. (2014a, b); Lawlor and Ferreira (2015) in order to calculate the forces acting between adatoms. Also non-interacting tight-binding model LeBohec et al. (2014), and Density Functional Theory (DFT) Solenov et al. (2013); Huang et al. (2012) were used in subsequent calculations. However, the influence of electron-electron interactions on the RKKY potential was not studied so far despite the fact that they are quite noticeable in graphene. Even the DFT approach is known to under-estimate the effect of inter-electron interactions. For instance, it strongly under-estimates the gap size in hydrogenated graphene Putz et al. (2014), which is strongly enhanced by interactions even at moderate concentrations of adatoms, as suggested by the recent QMC study in Ref. Ulybyshev and Katsnelson (2015). The importance of the effects of inter-electron interaction was also discussed in Black-Schaffer (2010).

In this paper we report on a first-principles Quantum Monte-Carlo (QMC) study of the RKKY interaction between hydrogen adatoms in graphene, consistently taking into account inter-electron interactions. We first consider pairwise interactions and demonstrate that for small distances between adatoms interaction effects dominate over the effects of finite temperature and finite hybridization terms. Then we consider the stability of adatom superlattices with respect to small shifts of a single adatom, finding the conditions for a dynamic stability of various superlattices with adatoms occupying only one or equally both sublattices. It will be shown that the stability conditions are substantially different from those previously defined in Ref. Abanin et al. (2010). The effect of antiferromagnetic (AFM) and charge density wave (CDW) ordering is also discussed, revealing an important feature of CDW order: the possibility to stabilize the superlattice configurations in which only one sublattice is occupied by hydrogen adatoms.

II Numerical setup.

We describe electrons in the conduction band of graphene using the standard tight-binding Hamiltonian with nearest-neighbor hoppings on the hexagonal lattice and electrostatic inter-electron interactions:

[TABLE]

where $\sum_{\langle x,y\rangle}$ and $\sum_{x,y}$ denote summations over all pairs $\langle x,y\rangle$ of nearest-neighbor sites and over all sites $x$, $y$ of the graphene honeycomb lattice respectively. $\hat{a}^{{\dagger}}_{x,\sigma}$, $\hat{a}_{x,\sigma}$ are the creation/annihilation operators for electrons with spin $\sigma=\uparrow,\downarrow$ in carbon $\pi$-orbitals, $t_{xy}$ are hopping amplitudes, $\hat{q}_{x}=-1+\sum_{\sigma}\hat{a}^{{\dagger}}_{x,\sigma}\hat{a}_{x,\sigma}$ is the charge operator at site $x$ and $V_{xy}$ is the inter-electron interaction potential. Some of the results presented in this work concern the non-interacting limit ($V_{xy}=0$), in case of which the model (1), with or without adatoms, can be solved exactly.

Periodic spatial boundary conditions are imposed as in Refs. Buividovich and Polikarpov (2012); Ulybyshev et al. (2013); Smith and von Smekal (2014):

[TABLE]

where $L_{1}$ and $L_{2}$ define the spatial size of the lattice.

In this work we use two models of hydrogen adatoms on the graphene sheet. The first is the simple vacancy model describing hydrogen adatoms as missing lattice sites in the tight-binding Hamiltonian (1), so that hopping amplitudes $t_{xy}$ are equal to zero for all neighbors $y$ of the lattice sites $x$ to which adatoms are attached. Away from adatoms, all hopping amplitudes are $t_{xy}=t=2.7\,{\rm eV}$. Furthermore, we assume that each adatom has zero charge.

The second model is the full hybridization model Wehling et al. (2010), in which hybridization terms

[TABLE]

are added to the tight-binding Hamiltonian (1), where $\gamma=2.0\,t$ is the hybridization parameter, $E_{d}=-0.06\,t$ is the electron energy for the adatom, $\sum_{x\in\mathcal{H}}$ denotes summation over all lattice sites with hydrogen adatoms and $\hat{c}^{{\dagger}}_{x,\sigma}$ are the creation operators for electrons on adatoms. As the hybridization model suffers from a fermion-sign problem (discussed in more detail below), which prevents the application of QMC, we use the full hybridization model only in the non-interacting limit. It is used, among other things, to verify the validity of the vacancy model and in cases where the effect of interactions can be modeled by an explicit anti-ferromagnetic or charge density wave (CDW) mass term.

In order to treat inter-electron interactions, we use the Suzuki-Trotter decomposition followed by the Hubbard-Stratonovich transformation and represent the partition function $\mathcal{Z}=\exp(-\hat{H}/T)$ at temperature $T$ as a path integral over Hubbard-Stratonovich fields $\phi_{x,\tau}$ in Euclidean time $\tau\in\left[0,T^{-1}\right]$ which is discretized, see Refs. Ulybyshev et al. (2013); Smith and von Smekal (2014) (to maintain exact particle-hole and sublattice symmetries one can use an exponential transfer matrix for the fermions between adjacent time slices Buividovich et al. (2016)):

[TABLE]

where $M_{e}=\partial_{\tau}-h_{xy}-i\phi_{x,\tau}\delta_{xy}$ and $M_{h}=\partial_{\tau}-h_{xy}+i\phi_{x,\tau}\delta_{xy}$ are the fermionic operators for electrons and holes respectively. The matrix of the single-particle tight-binding Hamiltonian $h_{xy}$ is identical for electrons and holes unless we introduce hybridization (4) or an additional mass term modelling a charge density wave. If the matrix $h_{xy}$ is the same, fermionic determinants for electrons and holes are complex conjugate:

[TABLE]

and the weight for the Hubbard fields in (5) is real and positive, which is a necessary requirement for a stochastic sampling of $\phi_{x,\tau}$. That this is not true in case of hybridization is the principle reason why we can use only the vacancy model in QMC calculations. Fortunately, as is demonstrated in the next section, this describes hydrogen adatoms with reasonable precision.

We sample the fields $\phi_{x,\tau}$ with the (manifestly positive) weight proportional to the integrand in (5) using the Hybrid Monte-Carlo algorithm. For the inter-electron interaction potential $V_{xy}$ we use the potentials calculated with the constrained RPA method Wehling et al. (2011) for suspended graphene (see Ulybyshev and Katsnelson (2015); Smith and von Smekal (2014) for details).

Within the interacting tight-binding model the RKKY interaction is nothing but the fermionic Casimir potential. For a pair of adatoms we calculate it as the free energy $\mathcal{F}_{xy}$ of the electrons on the graphene lattice with adatoms at sites $x$ and $y$ Bordag et al. (2001); Shytov et al. (2009). In absence of inter-electron interactions we simply compute the corresponding single-particle energy levels $\epsilon_{xy}$ with adatoms numerically and obtain $\mathcal{F}_{xy}$ up to an irrelevant constant $F_{0}$ from

[TABLE]

The free energy cannot be calculated directly in Hybrid Monte-Carlo simulations. To overcome this, we calculate the differences $\Delta\mathcal{F}=\mathcal{F}_{x+l,y}-\mathcal{F}_{x,y}$ between free energies for adatom positions which differ by a shift along one carbon-carbon lattice bond $l$. We represent this difference as an integral

[TABLE]

where $M_{\alpha}$ linearly interpolates between fermionic operators with adatoms at positions $x$ and $y$ (at $\alpha=0$) and $x+l$ and $y$ (at $\alpha=1$). Differentiating the path integral (9) for $\mathcal{Z}_{\alpha}$ by $\alpha$, we express $\Delta\mathcal{F}$ as

[TABLE]

where the expectation value is calculated with the same path integral weight as in (9). The matrix $\partial_{\alpha}M_{\alpha}$ is very sparse, allowing for an efficient calculation of ${\rm Tr}\,\left(M_{\alpha}^{-1}\partial_{\alpha}M_{\alpha}\right)$. The integral over $\alpha$ is calculated using the 6-point quadrature rule with six values of $\alpha\in\left[0,1\right]$, including $\alpha=0$ and $\alpha=1$. The above can be extended with no additional complications to cases with more than two adatoms.

III Pairwise interactions.

To study inter-adatom interactions in QMC simulations we use the simple vacancy model, since QMC has a fermion-sign problem for the hybridization model (4). For hydrogen adatoms the hybridization parameter $\gamma^{2}\gg E_{d}t$ is sufficiently large, so that the corresponding sp3 state of the carbon atom is effectively unavailable for $p_{z}$ electrons Katsnelson (2012) and the simple vacancy model is a good approximation to (4). In Fig. 1\textcolorreda we demonstrate that without inter-electron interactions the RKKY potentials are very similar for the hybridization model (4) and the vacancy model. In all cases, the pairwise interaction has well-known features: alternating signs for different sublattices Shytov et al. (2009) and the order-of-magnitude enhancement at some distances (clearly seen in Fig. 1\textcolorredb), at which the two adatoms induce midgap states with zero energy Solenov et al. (2013); LeBohec et al. (2014).

With reasonable computational resources, QMC simulations are limited to rather high temperatures in physical units ($T=0.09\,{\rm eV}=1040\,{\rm K}$ here) which are still relatively small, however, compared to the typical energy scales of the interacting tight-binding model in Eq. (1). In Fig. 1 we demonstrate that at least in absence of inter-electron interactions such temperatures indeed do not affect the qualitative features of the pairwise RKKY interactions. This suggests that we may safely use the vacancy model for adatoms in QMC at $T=0.09\,{\rm eV}=t/30$. We use lattices with $24\times 24$ cells, for which finite-volume effects are smaller than statistical errors.

In Fig. 2 we illustrate the effect of inter-electron interactions on the RKKY potential along the zigzag direction by comparison with the free electrons. In the free case we use the same vacancy model and exact numerical computation of the free energy according to eq. (7). The potential is particularly strongly modified at distances of 3 and 4 C-C bonds, while at distances larger then 8-9 C-C bonds the change of potential is too small to detect it with QMC. The main physical effect of the electron-electron interactions is that the local minimum at a distance of 3 bonds disappears and the potential barrier between widely separated adatoms and the global minimum corresponding to a dimer configuration becomes harder to penetrate.

IV Superlattices.

IV.1 Interaction effects

We now consider superlattices of regularly distributed hydrogen adatoms as examples of functionalization with a finite adatom concentration. To address superlattice stability, we consider the variation of free energy accompanying the shift of a single adatom from its regular position. First we compare interacting and non-interacting profile of the free energy accompanying the shift of one adatom along the zigzag direction. This profiles are shown on Figs. 3\textcolorreda and 3\textcolorredb both for interacting and non-interacting tight-binding models. We chose the system with 5.56% coverage of hydrogen adatoms populating only one sublattice, as illustrated on Fig. 3\textcolorredc. The vacancy model is used in both the interacting and non-interacting case (to avoid a sign problem for the former and make a direct comparison meaningful).

First we note that the overall scale of the RKKY interaction for superlattices is enhanced in comparison with pairwise interactions, so that the RKKY potential for a single adatom (with all other adatoms fixed) becomes comparable with the diffusion barriers $\Delta U\sim 0.3\ldots 1.0\,{\rm eV}$ for hydrogen adatoms Talbot et al. (2016).

Surprisingly, for superlattices inter-electron interactions do not change the RKKY potential qualitatively, despite inducing a very large gap $\Delta\epsilon\sim 1\,{\rm eV}$ in the midgap energy band Ulybyshev and Katsnelson (2015). To understand this observation, we recall that inter-electron interactions induce global anti-ferromagnetic (AFM) ordering for graphene with adatoms Ulybyshev and Katsnelson (2015), with the effective mass term $\hat{M}_{AFM}=\pm m\sum_{x}(\hat{a}^{{\dagger}}_{x\uparrow}\hat{a}_{x\uparrow}-\hat{a}^{{\dagger}}_{x\downarrow}\hat{a}_{x\downarrow})$ which has alternating signs on different sublattices. We can estimate the change in the free energy upon the shift of a single adatom to the neighboring lattice site (and thus to another sublattice) assuming that 1) this change is determined mostly by Tamm states localized near this adatom and 2) the energy of this Tamm state can be estimated as the mass term at nearest-neighbour sites of the adatom (as the wavefunctions of the Tamm states are mostly localized on these sites). Since the AFM mass has different signs for different spin components, the states with different spins simply exchange their energies, and the overall sum of energies in (7) doesn’t change in such a “mean field” approximation (see Fig. 4 for illustration).

The same argument applied to charge density wave (CDW) ordering leads to a completely different conclusion. Since in this case the mass term has the same sign for both spin components, the energies of two spin components no longer compensate each other, and the change of free energy upon adatom shift can be estimated as $2m$ (see Fig. 4). Unfortunately, verification of this scenario in full QMC simulations is very difficult, since a CDW mass term causes a sign problem: the fermionic determinant in (5) is no more positive definiteBuividovich et al. (2016). But at least we can illustrate the effects of both CDW and AFM mass terms on the RKKY potential for free electrons. The results are shown in Figs. 3\textcolorreda and 3\textcolorredb. We use $m=0.5\,{\rm eV}$, which approximately corresponds to the AFM mass induced by inter-electron interactions for this concentration of defects Ulybyshev and Katsnelson (2015). We indeed see that while the non-interacting result with AFM mass $m=0.5\,{\rm eV}$ almost coincides with the QMC result, the CDW mass term completely changes the RKKY potential and the locations of its minima.

We note, however, that the energy gap size of $0.5\,{\rm eV}$ is an over-estimate for the real graphene on boron nitride substrate. To demonstrate this, we have performed Monte-Carlo simulations on $6\times 6$ lattices with the bare CDW mass term $m_{CDW}=0.04$ eV which is close to that in real graphene on boron nitride Giovannetti et al. (2007). At nonzero $m_{CDW}$ the fermionic determinants for particles and holes are no longer complex conjugate to each other, thus the relation (6) is no longer valid and the path integral weight in the partition function (5) acquires some complex phase. We treat this complex phase using the brute-force reweighting. Namely, we sample the configurations of Hubbard-Stratonovich fields with the probability proportional to the absolute value of the path integral weight, and include the complex phase into the expectation values of physical observables like (10). Due to the fact that changing AFM mass term to the CDW mass term of the same value changes only the complex phase of the fermionic determinant but does not the absolute value, we were still able to represent the absolute value of the product of two determinants in (5) as the square of the absolute value of a single (electron or hole) determinant, and apply the same Hybrid Monte-Carlo algorithm as for the AFM mass term. We used the same setup as for the superlattice shown in the plots on the figures \textcolorred3a - \textcolorred3c, but reduced overall lattice size to $6\times 6$. Already at this small lattice size the complex phase exhibits strong oscillations which require a lot of statistics for reweighting. To illustrate these difficulties, on Fig. 5 we demonstrate that the distribution of the complex phase in the path integrals (5) is nearly flat, thus phase cancellations between different configurations are very important and require very large number of Monte-Carlo samples to resolve with good statistical accuracy. These difficulties in reweighting limit the simulations to rather small values of the CDW mass term $m_{CDW}$ and to small lattice sizes.

Nevertheless, even on small lattices we can obtain a rough estimate of the influence of CDW mass term induced by boron nitride on RKKY interaction potential. Using the reweighting technique we computed the change of the energy accompanying the shift of one adatom from it’s regular superlattice position to the nearest neighbor site. This shift corresponds to the position 1 in the plots on fig. \textcolorred3a and \textcolorred3b. The change of the energy is $\Delta F=-0.61\pm 0.16$ eV in the case of CDW mass while for zero mass the same calculation yields $\Delta F=-0.506\pm 0.018$ eV. We conclude that such a small bare CDW mass term is not enough to change the sign of the $\Delta F$. More generally, this means that real graphene is rather far from CDW phase transition so that the renormalization of the corresponding mass term due to inter-electron interactions is not very significant. On the other hand, the appearance of a large CDW mass term could still be expected if the ratio between on-site and nearest-neighbor electrostatic interaction potentials could be tuned to favor the CDW ordering Buividovich et al. (2016).

IV.2 Dynamic stability of superlattices

In order to address the dynamic stability of superlattice configurations with only one or both sublattices populated by adatoms, we use the hybridization model (4) without inter-electron interactions. This seems fairly well justified because these interactions do not qualitatively change the RKKY potential and their quantitative effect can be mimicked by an explicit AFM mass quite well (see Fig. 3\textcolorredb). Previously this subject was studied in the papers Abanin et al. (2010); Lawlor et al. (2014a, b); Lawlor and Ferreira (2015). However, the free energy of the system with defects was calculated in Lawlor et al. (2014a, b); Lawlor and Ferreira (2015) using some approximation for the fermionic propagator (Stationary Phase approximation) and for the free energy itself. For instance, only pairwise interactions were taken into account in the Monte Carlo study of the dynamic stability of various spatial configurations of defects in Lawlor and Ferreira (2015). Moreover, in all these papers, the calculation of the full free energy was used to study which spatial configuration of adatoms is energetically favourable. Randomly generated adatom configuration with equally/unequally populated sublattices were used in those studies. But these calculations don’t in general imply the stability of superlattices with respect to adatom displacements. The reason is that the real global energy minimum is collection of dimers due to very large pairwise RKKY interaction at nearest-neighbour position (see figures \textcolorred1 and \textcolorred2). So that there is no guarantee that a given spatial configuration of adatoms won’t change into a collection of dimers after set of energetically favourable shifts of adatoms’ positions. For this reason we study the change of the energy of superlattice after shift of one adatom, which automatically changes the relative occupation of sublattices.

Results of our calculations are presented in Fig. 6. We observe that the superlattices with adatoms on a single sublattice at half filling are dynamically unstable in all cases considered here due to fact that change of position of adatom to the opposite sublattice is energetically favourable in any case. In contrast, superlattices of adatoms which equally populate both sublattices are stable for low adatom concentration (see the structure in Fig. 6\textcolorredb). For higher adatom concentrations, all superlattices become unstable. This instability is likely to lead to the formation of dimers with a large binding energy.

Since a finite chemical potential can change the sublattice preferences of the pairwise RKKY interaction Solenov et al. (2013); LeBohec et al. (2014); Huang et al. (2012), the stability region of superlattices with only one sublattice occupied by adatoms might start from some finite chemical potential rather than at half-filling. In Fig. 7 we show the minimal change of the free energy among the three possible nearest-neighbor shifts as a function of the chemical potential $\mu$ for the same 6 adatom configurations as in Fig. 6. One can see that the single-sublattice superlattice with the lowest concentration of $3\%$ adatoms (label “a”) does indeed become stable above $\mu\approx 0.25\,{\rm eV}$. At about the same value of $\mu$ the corresponding superlattice with equal population of adatoms on both sublattices on the other hand becomes unstable. The structures with equally populated sublattices are stable mainly near half-filling for low concentrations of adatoms ($\leq 3\%$) or in some region of finite dopings for larger adatom concentrations. The larger the concentration, the smaller is the stability region. And vice versa for the single-sublattice configurations, the denser these adatom configurations get, the larger a chemical potential is needed to stabilize them.

Similar change can be observed in pairwise RKKY interactions, where positions of barriers and local minima interchange after some value of chemical potential. This phenomenon is illustrated in Fig. 8 for the same model of free electrons with hybridization term (4). The only difference that the critical chemical potentials seems to be much smaller (0.2 eV vs 0.4 eV) for superlattices.

V Conclusions.

We have calculated the RKKY interaction potential between hydrogen adatoms on a graphene sheet, taking into account effects of electron-electron interactions in fully non-perturbative first-principles QMC simulations. In particular, we have studied both, pairwise potentials and free-energy differences with stability analyses for various configurations of finite adatom densities. For the pairwise RKKY potential we found that the inter-electron interactions tend to increase the potetial barrier between widely separated adatom and dimer configurations which implies some suppression of dimer formation in the process of random deposition of adatoms an a graphene sheet.

For finite adatom concentrations, we have demonstrated that charge-density formation (CDW) and anti-ferromagnetic order (AFM) in the ground state, whether induced by substrates leading to staggered on-site potentials or dynamically by the inter-electron interactions, have very different effects. While an AFM mass term does not qualitatively change the RKKY-type interaction, the effect of a CDW mass term can be much more significant and even influence the sublattice ordering of adatoms.

Our stability analyses of different hydrogen superlattices show that single-sublattice configurations of adatoms are unstable at half filling but can be stabilized by an appropriate amount of doping with chemical potentials $|\mu|>\mu_{c}$. The critical value $\mu_{c}$ thereby increases with increasing adatom concentrations. Superlattices with equally populated sublattices are stable near half-filling for low concentrations of adatoms. More densely populated superlattices are likely to be unstable towards dimer formation.

As further plans we would like also to mention the variety of rich RKKY-related physics in bilayer graphene Klier et al. (2016). Taking into account that interaction effects might be also important in bilayer graphene Feldman et al. (2009), similar calculations for it are in our plans for future work.

Acknowledgements.

We thank M. I. Katsnelson for useful and motivating discussions. This work was supported by the Deutsche Forschungsgemeinschaft (DFG), grants BU 2626/2-1 and SM 70/3-1. Calculations have been performed on GPU clusters at the Universities of Giessen and Regensburg. P.B. was supported by the S. Kowalevskaja award from the Alexander von Humboldt foundation.