{'fraction_non_alphanumeric': 0.03119498915134892, 'fraction_numerical': 0.010275514799197609, 'mean_word_length': 4.849616858237548, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Inter-electron interactions and the RKKY potential between H adatoms in graphene

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.


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 [1][2][3] with a possibility to tune it and even to return the material to the initial semi-metallic state [4]. Also the magnetic moments induced around hydrogen adatoms due to inter-electron interactions [3,5,6] play an important role in spin relaxation processes [7] 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 [3]. 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 [8,9]. Especially important is the case of functionalized graphene on top of boron nitride [8], 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 ∼ 10 meV for distances larger than the inter-atomic spacing [10] and decays as r −3 at large distances [11], suggests that the Ruderman-Kittel-Kasuya-Yosida (RKKY) con-tribution 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 [12]. Typically, some analytic approximations to the fermionic Green's function in presence of resonant scatterers are used [13][14][15][16] in order to calculate the forces acting between adatoms. Also non-interacting tight-binding model [17], and Density Functional Theory (DFT) [18,19] 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 [2], which is strongly enhanced by interactions even at moderate concentrations of adatoms, as suggested by the recent QMC study in Ref. [3]. The importance of the effects of inter-electron interaction was also discussed in [20].
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. [1]. 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.

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: where x,y and x,y denote summations over all pairs x, y of nearest-neighbor sites and over all sites x, y of the graphene honeycomb lattice respectively.â † x,σ ,â x,σ are the creation/annihilation operators for electrons with spin σ =↑, ↓ in carbon π-orbitals, t xy are hopping amplitudes,q x = −1 + σâ † x,σâx,σ 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. [21][22][23]: 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 eV. Furthermore, we assume that each adatom has zero charge.
The second model is the full hybridization model [24], in which hybridization termŝ are added to the tight-binding Hamiltonian (1), where γ = 2.0 t is the hybridization parameter, E d = −0.06 t is the electron energy for the adatom, x∈H denotes summation over all lattice sites with hydrogen adatoms and c † x,σ 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 Z = exp(−Ĥ/T ) at temperature T as a path integral over Hubbard-Stratonovich fields φ x,τ in Euclidean time τ ∈ 0, T −1 which is discretized, see Refs. [22,23] (to maintain exact particle-hole and sublattice symmetries one can use an exponential transfer matrix for the fermions between adjacent time slices [25]): where M e = ∂ τ − h xy − iφ x,τ δ xy and M h = ∂ τ − h xy + iφ x,τ δ xy are the fermionic operators for electrons and holes respectively. The matrix of the single-particle tightbinding 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: (6) and the weight for the Hubbard fields in (5) is real and positive, which is a necessary requirement for a stochastic sampling of φ x,τ . 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 φ x,τ 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 [26] for suspended graphene (see [3,23] 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 F xy of the electrons on the graphene lattice with adatoms at sites x and y [13,27]. In absence of interelectron interactions we simply compute the corresponding single-particle energy levels ǫ xy with adatoms numerically and obtain F xy up to an irrelevant constant F 0 from F xy = −T ǫxy ln 1 + e −ǫxy/T + F 0 .
The free energy cannot be calculated directly in Hybrid Monte-Carlo simulations. To overcome this, we calculate the differences ∆F = F x+l,y − 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 where M α linearly interpolates between fermionic operators with adatoms at positions x and y (at α = 0) and x + l and y (at α = 1). Differentiating the path integral (9) for Z α by α, we express ∆F as where the expectation value is calculated with the same path integral weight as in (9). The matrix ∂ α M α is very sparse, allowing for an efficient calculation of The integral over α is calculated using the 6-point quadrature rule with six values of α ∈ [0, 1], including α = 0 and α = 1. The above can be extended with no additional complications to cases with more than two adatoms.

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 γ 2 ≫ E d t is sufficiently large, so that the corresponding sp 3 state of the carbon atom is effectively unavailable for p z electrons [28] and the simple vacancy model is a good approximation to (4). In Fig. 1a 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 [13] and the order-of-magnitude enhancement at some distances (clearly seen in Fig. 1b), at which the two adatoms induce midgap states with zero energy [17,18].
With reasonable computational resources, QMC simulations are limited to rather high temperatures in physical units (T = 0.09 eV = 1040 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 interelectron 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 eV = t/30. We use lattices with 24 × 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. 4. SUPERLATTICES.

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 noninteracting profile of the free energy accompanying the shift of one adatom along the zigzag direction. This profiles are shown on Figs. 3a and 3b 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. 3c. The vacancy model is used in both the interacting and noninteracting 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 ∆U ∼ 0.3 . . . 1.0 eV for hydrogen adatoms [29].
Surprisingly, for superlattices inter-electron interactions do not change the RKKY potential qualitatively, despite inducing a very large gap ∆ǫ ∼ 1 eV in the midgap energy band [3]. To understand this observation, we recall that inter-electron interactions induce global anti-ferromagnetic (AFM) ordering for graphene with adatoms [3], with the effective mass termM AF M = ±m x (â † x↑â x↑ −â † x↓â x↓ ) 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 com-  ponents 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 definite [25]. 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. 3a and 3b. We use m = 0.5 eV, which approximately corresponds to the AFM mass induced by inter-electron interactions for this concentration of defects [3]. We indeed see that while the non-interacting result with AFM mass m = 0.5 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 eV is an over-estimate for the real graphene on boron nitride substrate. To demonstrate this, we have performed Monte-Carlo simulations on 6 × 6 lattices with the bare CDW mass term m CDW = 0.04 eV which is close to that in real graphene on boron nitride [30]. 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 3a -3c, but reduced overall lattice size to 6 × 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. 3a and 3b. The change of the energy is ∆F = −0.61 ± 0.16 eV in the case of CDW mass while for zero mass the same calculation yields ∆F = −0.506 ± 0.018 eV. We conclude that such a small bare CDW mass term is not enough to change the sign of the ∆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 [25].

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. 3b). Previously this subject was studied in the papers [1,[14][15][16]. However, the free energy of the system with defects was calculated in [14][15][16] 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 [16]. 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 1 and 2). 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. 6b). 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 [17][18][19], 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 µ 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 µ ≈ 0.25 eV. At about the same value of µ 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 (≤ 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.

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 interelectron 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 antiferromagnetic 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 |µ| > µ c . The critical value µ 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 [31]. Taking into account that interaction effects might be also important in bilayer graphene [32], similar calculations for it are in our plans for future work.