A theoretical study of the dynamics of atomic hydrogen on graphene bilayers

A theoretical study of the dynamics of atomic hydrogen on graphene bilayers

Mohammed Moaied moaied5@yahoo.com Departamento de Física de la Materia Condensada, Universidad Autónoma de Madrid, Cantoblanco, 28049 Madrid, Spain and Department of Physics, Faculty of Science, Zagazig University, 44519 Zagazig, Egypt.    J. A. Moreno joseantonio.moreno@estudiante.uam.es Departamento de Física de la Materia Condensada, Universidad Autónoma de Madrid, Cantoblanco, 28049 Madrid, Spain.    M. J. Caturla mj.caturla@ua.es Departamento de Física Aplicada, Universidad de Alicante, San Vicente del Raspeig, 03690 Alicante, Spain.    Félix Ynduráin    J. J. Palacios juanjose.palacios@uam.es Departamento de Física de la Materia Condensada, Instituto Nicolás Cabrera (INC), and Condensed Matter Physics Institute (IFIMAC), Universidad Autónoma de Madrid, Cantoblanco, 28049 Madrid, Spain.
July 16, 2019

We present a theoretical study of the dynamics of H atoms adsorbed on graphene bilayers with Bernal stacking. First, through extensive density functional theory calculations, including van der Waals interactions, we obtain the activation barriers involved in the desorption and migration processes of a single H atom. These barriers, along with attempt rates and the energetics of H pairs, are used as input parameters in kinetic Monte Carlo simulations to study the time evolution of an initial random distribution of adsorbed H atoms. The simulations reveal that, at room temperature, H atoms occupy only one sublattice before they completely desorb or form clusters. This sublattice selectivity in the distribution of H atoms may last for sufficiently long periods of time upon lowering the temperature down to C. The final fate of the H atoms, namely, desorption or cluster formation, depends on the actual relative values of the activation barriers which can be tuned by doping. In some cases a sublattice selectivity can be obtained for periods of time experimentally relevant even at room temperature. This result shows the possibility for observation and applications of the ferromagnetic state associated with such distribution.

I Introduction

Hydrogenation of carbon-based materials such as graphene, graphite, or carbon nanotubes is attracting much interest as a practical methodology to manipulate the electronic and magnetic properties of these materials in a reversible way. Hydrogenation of graphene, for instance, was found, both theoretically and experimentally, to be an effective way to turn this system from a gapless semiconductor into a gapful one with a tunable band gap(Elias et al., 2009; Sessi et al., 2009; Haberer et al., 2010; Yang et al., 2010; Balog et al., 2010; Lin et al., 2015). Controlled hydrogenation, on the other hand, has been predicted to induce interesting magnetic states with potential applications in spintronics (Zhou et al., 2009; Soriano et al., 2010; Leconte et al., 2011; McCreary et al., 2012; Soriano et al., 2011; Nair et al., 2012). Unfortunately, the necessary control has not been experimentally demonstrated to date.

Graphene is a single layer of carbon atoms bonded together in a bipartite honeycomb structure which is formed by two inter-penetrating triangular sublattices. The first neighbors of an atom in a given sublattice belong to the other sublattice and vice versa (Saito et al., 1998). When atomic H is adsorbed on graphene the H atom bonds directly on top of a carbon atom and induces an intrinsic magnetic moment around the adsorption site with a net magnetic moment of 1 (Yazyev and Helm, 2007; Boukhvalov et al., 2008; Casolo et al., 2009; Soriano et al., 2010; Yndurain, 2014). Since the sublattices are chemically equivalent, the adsorption process is blind to the sublattice index. For graphite or multilayer graphene, unlike the monolayer case, surface carbon atoms in one sublattice present others underneath while the ones in the complementary sublattice do not (assuming Bernal stacking). Thus, at first sight, it should not come as a surprise that the two sublattices offer different binding energies to H atoms, thus favouring adsorption on one of them (at least at low enough coverages).

The specific details of the desorption and diffusion processes of the adsorbed H atoms are, nevertheless, essential to determine the final or temporary hydrogenation patterns and related electronic properties and, most importantly, the time scale for reaching such patterns. When H atoms are initially deposited, e.g., by cracking molecular HGüttler et al. (2004), it is expected that they will reach both sublattices with equal probability. The electronic state thus induced will correspond to that of a non-magnetic insulator for large concentrations(Elias et al., 2009; Sessi et al., 2009; Haberer et al., 2010; Yang et al., 2010; Balog et al., 2010), an antiferromagnet for intermediate concentrationsPalacios et al. (2008), or a paramagnet for low concentrationsNair et al. (2012). At room temperature, which for practical applications is the most interesting case, both desorption and diffusion processes are, in principle, active(Jeloaica and Sidis, 1999; Hornekær et al., 2006; Chen et al., 2007; Denis and Iribarne, 2009; Kerwin and Jackson, 2008; Ivanovskaya et al., 2010; Sha and Jackson, 2002; Herrero and Ramirez, 2009; Huang et al., 2011a, b; Borodin et al., 2011). If desorption rates are larger than diffusion ones, the sample will loose H from both sublattices, but at different rates. If, on the other hand, diffusion or migration rates are larger than desorption ones, H atoms will move across the surface performing many jumps before desorbing. In this case they will certainly spend more time on one sublattice than on the other. In both scenarios one sublattice may become more populated than the other, although in the second one the H atoms will also have more chances to come across one another and form stable non-magnetic clusters(Ferro et al., 2003; Casolo et al., 2009; Sljivancanin et al., 2009; Roman et al., 2007; Andree et al., 2006).

Here we show that one can get, at least temporarily, a nearly 100% sublattice selectivity in the adsorption, i.e., a distribution where all adsorbed H atoms occupy the same sublattice on the surface graphene monolayer. At room temperature this may occur for periods of time of minutes. Upon lowering the temperature down to C, the single-sublattice distribution may, however, survive for days. This result does not qualitatively depend on the specific values of the migration and desorption barriers. For instance, upon changing the carrier concentration of the bilayer system, which, in turn, changes the magnitude of these barriersHuang et al. (2011b), we always obtain such temporary distribution of H, only the final fate of the atoms being affected. For hole doping, all atoms eventually form dimers or clusters while for electron doping (or no doping) all atoms eventually desorb. This remarkable result has important implications since theory predicts that when all the H atoms bind to the same sublattice the resulting electronic ground state should be a ferromagnetic state(Lieb, 1989; Zhou et al., 2009; Palacios et al., 2008) with a typically very high Curie temperature for a wide range of concentrationsMoaied et al. (2014).

Ii Computational Methodology

We combine two types of computational methodologies. On the one hand, density functional theory (DFT)(Hohenberg and Kohn, 1964; Kohn and Sham, 1965) and, on the other, kinetic Monte Carlo (KMC) simulations. The essential ingredients for the latter, i.e., activation energies and attempt rates, are obtained from the former. A graphene bilayer is considered all along; the results thus obtained need not be exactly representative of the physics of multilayer graphene or graphite, although arguments can be put forward to carry our results over to these cases.

ii.1 van der Waals DFT with SIESTA and preliminary checks

As we are dealing with weakly interacting graphene layers where dispersion (van der Waals) forces due to long-range electron correlation effects play a key role, we employ the non-local van der Waals density functional (vdW-DF) of Dion et al. (Dion et al., 2004) as implemented by Román-Pérez and Soler(Román-Pérez and Soler, 2009; Kong et al., 2009) in the SIESTA code(Ordejon et al., 1996; Soler et al., 2002). To describe the interaction between the valence and core electrons we used norm-conserved Troullier-Martins pseudopotentials (Troullier and Martins, 1991). To expand the wave functions of the valence electrons, a double- plus polarization (DZP) basis set was used (Junquera et al., 2001). We experimented with a variety of LCAO basis sets and found that including polarization basis elements was important; however, a triple- (TZP) basis set produced results essentially indistinguishable from those obtained with DZP. For the Brillouin zone sampling we used at least a Monkhorst-Pack mesh, increasing the density of points for occasional checks. We ensured that the vacuum space is at least 25 Å so that the interaction between functionalized layers and their periodic images can be safely ignored. We have also checked that the results are well converged with respect to the real space grid. The atoms are allowed to relax down to a force tolerance of 0.020 eV/Å, keeping the necessary coordinates of the H atom fixed to obtain the corresponding desorption energy curves and migration landscapes. Spin polarization was included in the calculations because hydrogenation induces magnetism in single-layer as well as multilayer grapheneMoaied et al. (2014).

Figure 1: (color online). Binding energy of a H atom on top a C atom in monolayer graphene (a supercell) as a function of the distance to the graphene plane. Blue dots correspond to a standard GGA functional while black ones correspond to the vdW-DF, as explained in the text.

Before addressing the specifics of bilayer graphene in Sec. III, we will briefly revisit the monolayer graphene case. H atoms are known to preferentially adsorb on top of carbon atoms. Figure 1 shows the vdW-DF binding energy (black dots) of a H atom on top a C atom for a single graphene layer as a function of the distance to the graphene plane, . The binding energy is defined as usual:


(As a complementary accuracy check, we have made sure that the energy of the isolated H atom, , is 1 Ryd.) Two distinctive minima or adsorption states appear: a strongly bound chemisorption state at Å  with binding energy and a weakly bound physisorption state, which can be appreciated as a shallow minimum around Å. The distance between the host C atom and the H atom abruptly changes in between both minima, being Å  at the chemisorption state and Å  at the physisorption state. The derivative discontinuity at the transition point between energy minima can be attributed to the mean field treatment which could be smoothed out by more sophisticated methods which do not break spin symmetryCasolo et al. (2009).

The importance of using a vdW-DF not only reflects on the fact that standard GGA functionals do not bind (or barely bind) graphene layers into a bilayer or graphite. Blue dots in Fig. 1 correspond to the same calculation using a commonly used GGA functionalZhang and Yang (1998). The result is essentially similar to many others found in the literature(Jeloaica and Sidis, 1999; Hornekær et al., 2006; Chen et al., 2007; Denis and Iribarne, 2009; Kerwin and Jackson, 2008; Ivanovskaya et al., 2010; Sha and Jackson, 2002; Herrero and Ramirez, 2009; Huang et al., 2011a, b; Borodin et al., 2011), whereas it is quantitatively and even qualitatively different from the vdW-DF result. In particular, deeper chemisorption and physisorption minima are obtained in the vdW-DF calculation. Although the saddle point separating the chemisorption state from the physisorption one is not smooth in our numerics, we can still appreciate that it has a negative value meV. Most calculations in the literature exhibit slightly positive saddle point energies(Jeloaica and Sidis, 1999; Hornekær et al., 2006; Chen et al., 2007; Denis and Iribarne, 2009; Kerwin and Jackson, 2008; Ivanovskaya et al., 2010; Sha and Jackson, 2002). Interestingly, the difference is similar for both functionals.

ii.2 Object kinetic Monte Carlo algorithm

The KMC calculations have been performed to understand the time evolution of the surface distribution of the H atoms for different temperatures and concentrations. KMC algorithms are powerful techniques to study the dynamics of a system of particles when the different events that those particles can perform are known as well as their probabilities (for a recent review see (Voter, 2007)). There are many different algorithms with the name KMC. In this case we use what is often known as an object kinetic Monte Carlo (OKMC) algorithm, based on the residence time algorithm or Bortz-Kalos-Liebowitz (BKL) algorithm (Kalos and Whitlock, 1986; Voter, 2007) Briefly, in an OKMC algorithm a list of possible events is defined with a given probability for each event, . This probability usually follows an Arrhenius dependence with temperature:


where is the Boltzmann constant, is the activation energy of the given event, and is the attempt frequency. In this case the activation energies are related to migration or diffusion energies, desorption, and dissociation energies.

The total rate for all events, , is then calculated as:


where is the total number of events and is the number of particles that can perform event . An event is then selected by picking a random number between . In this way one event is selected every Monte Carlo step from all possible with the appropriate weight. Once the event has been selected, a random particle is chosen from all those that can undergo that event. The particle is then moved and the total rate has to be computed again for the next simulation step. At every Monte Carlo step the time increases by:


where is a random number between that is used to give a Poisson distribution of the time.

Figure 2: (color online). Atomic structure of H on bilayer graphene. (a) and (b) sites top view and (c) and (d) sites side view detail.

Iii Activation barriers for a H atom on bilayer graphene

iii.1 Desorption barriers

The physics of atomic H adsorption on a bilayer is not much different from that on monolayer graphene. Now only the upper layer is allowed to relax while the C atoms in the lower layer are being fixed at their equilibrium positions, simulating the presence of a substrate such as graphite or SiC. Due to the chosen Bernal stacking, the two sublattices are not equivalent any more. To stress this important point, we will denote the two different adsorption sites as and from now on (see Fig. 2). The vdW-DF binding energy curves corresponding to both adsorption sites are presented in Fig. 3 for a supercell. As for the monolayer case (see Fig. 1), both curves exhibit two minima or adsorption states. A magnetic moment of 1 also appears on the surface layer at the chemisorption state while it transfers to the H atom at the physisorption minimum. In Fig. 2 we show the atomic structure for the chemisorption state. The characteristic re-hybridization induced by the H atom is patent in both adsorption sites. The resulting atomic structures are not identical (although this can barely be appreciated in the figure) and the chemisorption energies, and , are slightly different as well (). No significant difference between adsorption sites can be appreciated in the physisorption part of the curves (see Fig. 3).

Similarly to the monolayer case, the saddle points separating chemisorption from physisorption minima, (), present negative values ( meV). In the light of this result one might wonder whether the desorption activation barriers, , to be considered in the OKMC calculations, should correspond to or to the smaller difference . Figure 4 shows a blow-up of the physisorption minima for the H atom being on top a site, hollow, and bridge positions. As can be appreciated, the physisorption energy barely changes with the position of the H atom with respect to the substrate; in other words, the physisorption energy landscape is essentially flat on the scale of eV. H atoms being thermally promoted from the chemisorption to the physisorption state across the saddle point may freely wander across the surface before they are chemisorbed again or finally desorb. In what follows we will assume that the distance travelled on the physisorption channel is long enough so that, before the chemisorption process occurs again, the H atoms meet other H atoms or hit grain boundaries where they permanently stick (see below). Thus, to any practical purpose, we will consider . (The possibility of considering as the relevant desorption activation barrier is briefly discussed below).

Figure 3: (color online). Binding energy curves for a H atom adsorbed on the (red) and (blue) sites of a bilayer graphene surface (a supercell).
Figure 4: (color online). Detail of the physisorption binding energy curve for a H atom on top a C atom ( site), on hollow, and on bridge position for a bilayer graphene (a supercell).

iii.2 Migration barriers

Figure 5 shows the binding energy obtained by displacing the H atom along the bond line () joining the and sites (for a supercell). These results are obtained by fixing the coordinates (plane) of the H atom, letting the coordinate (height) and all the positions of the carbon atoms of the top graphene layer to relax. It is important to notice that the sets of points originating from both sites (differentiated by colors) do not cross in the coordinates phase space. At the “crossing” point around Å, two different solutions with very different atomic structures are obtained, one where the H atom remains bonded to the C atom and the other one where the H atom is bonded to the C atom (see insets). In fact, as shown in the figure, both sets of points can be smoothly continued beyond the crossing point. While the actual path in coordinates phase space for the H atom to move from the site to the one (or viceversa) is not known to us, it must cross a saddle point where the distance between the H atom and the two C atoms involved is the same. The binding energy of this saddle point is meV (represented by a black dot in Fig. 5) and its atomic structure is shown in the upper inset. Taken with respect to the and chemisorption minima, the saddle point energy gives the activation barriers to directly migrate between sites, . From the figure we see that these barriers are essentially equal to the respective chemisorption energies , but are larger than the desorption barriers as defined above (by meV). Importantly, as shown in Fig. 4, the H atom is still considerably bound to the surface at the bridge or saddle point position, being thus possible to directly migrate between sites without moving to the physisorption channel or desorbing at room temperature.

Figure 5: (color online). Binding energy of a H atom on bilayer graphene placed along the bond that links the site to the site for a 4x4 supercell. Despite the appearances, the crossing point is not a saddle point (see lower insets and see text). The true saddle point is depicted with a black dot and the associated structure is shown in the upper inset. (In all insets the bottom layer is not shown for clarity).

iii.3 Activation processes involving two or more H atoms

The activation energies discussed above dramatically change when the H atom is close to other H atoms. Here we will make use of the energetics of two H atoms on a graphene monolayer (as reported, e.g., in Ref. Bachellerie et al., 2007; Rougeau et al., 2006; Sljivancanin et al., 2009; Moaied et al., 2014) to estimate how the vicinity of other atoms modify these barriers. A thorough study for some cluster possibilities can be found in Ref. Sljivancanin et al., 2009 which also serves of guidance for the following considerations. Our basic assumption is that when a H atom attempts to break an “bond” [two H atoms sitting on nearest-neighbour C atoms, also known as orto (O) dimerSljivancanin et al. (2009)] an extra pair binding energy eV has to be paidMoaied et al. (2014) in addition to the migration barrier calculated for isolated atoms. In other words, the energy required to change an H dimer into an or dimer [metastable (M) dimersSljivancanin et al. (2009)] is given by


The same assumption will be made for migration processes involving the breaking of pairs [also known as para (P) dimers (see Ref. Sljivancanin et al., 2009)], which are also strongly bound with a binding energy of eVMoaied et al. (2014). Finally, we also apply the same addition rule to the desorption barrier of a H atom if, in the desorption process, an or bond is broken:


In order not to complicate in excess and unnecessarily the OKMC simulations, additional assumptions regarding the energetics need to be made: a) We will ignore the binding energy or attraction between H atoms exerted at distances longer than those in the above referred pairs; b) H clusters (more than two atoms in close proximity) can always be considered as composed of dimers linked by and/or bonds so that when a H atom attempts a migration or desorption, the activation barrier will be that of breaking the corresponding bond(s), regardless of the cluster structure and number of H atoms forming it; c) finally, we will assume that the attempt rates are not modified by the presence of other H atoms. The accuracy of all these estimates is in fact not critical at all to the final results. Once clusters are formed, the activation energies are so large that they never break apart (in a relevant time scale and for relevant temperatures). When interpreting the results, one should only keep in mind that ignoring the interaction between H atoms at longer distances might reduce the likelihood for formation of clusters.

Event barrier frequency
Isolated atom events
Migration from site 1.15 3 3.5
Migration from site 1.23 3 3.5
Desorption of site 1.00 7.10
Desorption of site 1.08 7.10
Dissociative events
Migration from in dimer 2.55 2 3.5
Migration from in dimer 2.63 2 3.5
Migration from in dimer 2.5 3 3.5
Migration from in dimer 2.58 3 3.5
Desorption of in dimer 2.40 7.10
Desorption of in dimer 2.48 7.11
Desorption of in dimer 2.35 7.10
Desorption of in dimer 2.43 7.11
Table 1: Events included in the OKMC calculation along with the correspoding activation barriers (in eV) and attempt frequencies (in s) for each type of event, as obtained from the DFT calculations. The factor in front of the frequency values relates to the number of available positions to jump to.

Iv Kinetic Monte Carlo simulations

From simple statistical considerations, the fact that and chemisorption energies are different and that both migration and desorption activation barriers are in the vicinity of eV anticipate that near room temperature H atoms may occupy and sites with a significantly different probability. The question that actually needs to be addressed is how the two activated migration and desorption processes compete to determine the evolution of an initial random arrangement of H atoms and whether or not, in the end, all H atoms desorb or group together forming non-magnetic clusters, which, according to Eq. 5 should be thermodynamically stableŠljivančanin et al. (2012).

These questions can be answered through OKMC calculations. We have thus implemented an OKMC algorithm that includes a total of 12 events, as shown in Table 1. The rates for migration and desorption are estimated from the DFT activation energies, as explained above, and from the attempt frequencies, calculated using simple harmonic models derived from the energy curves near the minima. Despite of our efforts to obtain accurate activation energies, it is unrealistic to expect a precision down to meV or even tens of meV, so the overall time scales may still be somewhat uncertain. However, the difference between the and desorption and migration barriers is undeniable. This difference has been set to meV, which approximately corresponds to the chemisorption energy difference in the zero concentration limitMoaied et al. (2014). The results are thus expected to be generically representative of low concentrations. In all the simulations we have considered a number of lattice sites in the range of in order to have access to low values of the coverages (defined as the ratio of H atoms to C atoms) for a sufficiently large number of H atoms (). This keeps the statistical noise in the curves sufficiently low.

The actual chemisorption process after cracking molecular H is unknown to us and probably needs a separate discussion. In what follows we will assume an initial random distribution of chemisorbed H atoms for the chosen concentration. A H atom can either jump to a neighbouring location or desorb from the surface layer. After every jump it is necessary to check whether or not other H atoms are located in the vicinity. As explained above, we consider the formation of two types of pairs or dimers and clusters formed out of them which are kept immobile as a whole in the calculation. As a result we will show the time evolution of the relative abundance of H atoms adsorbed on each sublattice [ (in red) and (in blue)]. Solid lines will correspond to isolated atoms while H atoms belonging to dimers will be represented by dashed lines and those associated with dimers by dotted lines. In fact, since clusters with more than two H atoms can be formed, we should generally speak about H atoms being part of bonds instead of dimers. For instance, a trimer such as contains two atoms and one atom which we associate with bonds in our analysis. Therefore, as can be seen in some cases below, the number of and atoms associated with these bonds does not need to be identical. Likewise, the count of H atoms associated with dimers must be interpreted similarly. Cases such as atoms in clusters containing different types of bonds, e.g., are associated with both, but these seldom appear and do not significantly affect the count. We will also show the number of desorbed atoms (black lines).

Figure 6: (color online). Time evolution of the relative abundance of H atoms adsorbed on the two graphene sublattices for an initial concentration and two representative temperatures (a) T=300 K and (b) T=273 K. Solid lines correspond to isolated H atoms while dashed and dotted lines correspond to H atoms forming part of dimers (or clusters) through short or long bonds, respectively (see text for a detailed explanation of these two types of bonds). Black solid line refers to desorbed atoms.
Figure 7: (color online). Same as in Fig. 6, but only at room temperature for (a) and (b) .

iv.1 Neutral graphene bilayer

Representative results for different initial values of the H concentration and different temperatures for a neutral bilayer are presented in Figs. 6 and 7. In all cases we generically observe that the population of isolated H atoms on both sublattices decreases overtime due to desorption. Expectedly, after a certain time, 100% of the remaining H atoms sit on the sublattice due to the larger desorption barrier of this site. At room temperature it takes 1 hour for this to occur [see Fig. 6(a)] and after a few hours most of H has desorbed. At lower temperatures ( K) the time window to have a full concentration of H atoms on the same sublattice logically increases [to approximately several days as shown in Fig. 7(b)]. Interestingly, this single-sublattice distribution now lasts for hours before desorption takes over. In both cases cluster formation barely occurs because of the significant difference in the desorption and migration activation barriers. (However, we should keep in mind that we have excluded the possibility of migration on the physisorption channel which might increase the probability of cluster formation.) The whole picture remains essentially the same for a wide range of H concentrations, as shown in Fig. 7. Apart from an obvious statistical smoothing of the curves for larger concentrations, the probability for an initial accidental formation of dimers or clusters is, as expected, larger for larger , although the number of clusters does not change overtime.

Figure 8: (color online). Same as Fig. 6, but for a hole-doped bilayer. The carrier concentration has been chosen as to make equal the desorption and migration barriers. The H concentration is and the temperature is (a) 300 K and (b) 273 K.

iv.2 Doped graphene bilayer

We now examine how the dynamics is affected upon doping or carrier concentration variations. DFT calculations for migration and desorption activation barriers in a doped graphene monolayer have been reported in Ref. Huang et al., 2011b. There it is shown that doping affects both migration and desorption activation barriers. Following the results in Ref. Huang et al., 2011b and assuming that the doping reaches the upper layer if induced by the substrate or by a field-effect configuration, we have considered the case of hole doping where the desorption barrier becomes higher than the migration one (this already happens for very small concentration of carriersHuang et al. (2011b)). The dynamics in the case of electron doping is essentially similar to the one in the neutral case, only changing the time scales and will not be explored here.

In Fig. 8 we plot the time evolution of the relative H abundance for a hole concentration of cm at room temperature (a) and K (b). According to Ref. Huang et al., 2011b, this small doping changes the activation barriers so as to make the desorption and migration barriers alikeHuang et al. (2011b) in our case. Similarly to the neutral bilayer, after a certain time (which depends on temperature), all H atoms are hosted by the same sublattice. Now, however, migration starts playing a significant role since the desorption is hindered and the concentration of H atoms on the sublattice increases at the expense of the atoms initially adsorbed on the sublattice. Importantly, instead of desorbing at long times, some H atoms eventually form dimers or clusters (% for the chosen concentration of H). Incidentally, as discussed in Sec. III, had we considered the desorption barriers to be , the dynamics of H on a neutral bilayer would have been qualitatively similar to the one shown in Fig. 8.

In Fig. 9 we show the time evolution of the H abundance for a larger carrier concentration of cm, again at room temperature (a) and at a lower temperature (b). A smaller percentage of H atoms are lost now, the rest of them remaining adsorbed on the surface in the form of dimers or clusters at long times. On changing the initial concentration of H atoms (not shown) the ratio of desorbed to adsorbed atoms logically changes (the larger the concentration, the smaller the percentage of desorbed atoms). Interestingly, even at room temperature, there is now a time window from approximately 1 min to several days where a significant fraction of H atoms remain adsorbed on the same sublattice, coexisting with dimers or clusters at longer times. Notice that dimers are more abundant than ones, which can be understood since their formation probability is larger due to the larger capture radius. As in the previous cases, all time scales increase on lowering the temperature [see Fig. 9(b)].

Figure 9: (color online). Same as Fig. 8, but for a larger carrier concentration which increases even further the desorption barriers and decreases the migration barriers. Two different temperatures are shown: (a) K and (b) K. The concentration of H has been here set to .

V Conclusions and final considerations

Our combination of DFT calculations and OKMC simulations have shown that a selective sublattice adsorption of atomic H can be realized on the surface of graphene bilayers for low concentrations. Although, over time, H atoms eventually desorb and/or form clusters, depending on the doping conditions of the substrate, one can always find a time window where the selectivity persist and allows for measurements of the expected ferromagnetic properties of this system. The time window becomes large enough for measurements or even practical applications upon lowering the temperature just a few degrees down to C.

All the calculations presented here have been carried out for bilayer graphene. However, given the weak interaction between layers, no qualitative changes are expected for multilayer graphene or graphite. Finally, we should stress that, although the accuracy of DFT calculations is questionable down to the meV range, the uncertainty only affects overall times scales which, as shown, can be easily tuned with temperature or doping.

Vi Acknowledgments

This work was supported by MINECO under Grants Nos. FIS2013-47328 and FIS2012-37549, by CAM under Grants Nos. S2013/MIT-3007, P2013/MIT-2850, and by Generalitat Valenciana under Grant PROMETEO/2012/011. The authors acknowledge I. Brihuega for many discussions and for sharing unpublished experimental results with us and J. Soler for enlightening discussions. The authors thankfully acknowledge the computer resources, technical expertise, and assistance provided by the Centro de Computación Científica of the Universidad Autónoma de Madrid and by the Supercomputing and Visualization Center of Madrid (CeSViMa).


  • Elias et al. (2009) D. C. Elias, R. R. Nair, T. M. G. Mohiuddin, S. V. Morozov, P. Blake, M. P. Halsall, A. C. Ferrari, D. W. Boukhvalov, M. I. Katsnelson, A. K. Geim, et al., Science 323, 610 (2009).
  • Sessi et al. (2009) P. Sessi, J. R. Guest, M. Bode, and N. P. Guisinger, Nano letters 9, 4343 (2009).
  • Haberer et al. (2010) D. Haberer, D. V. Vyalikh, S. Taioli, B. Dora, M. Farjam, J. Fink, D. Marchenko, T. Pichler, K. Ziegler, S. Simonucci, et al., Nano letters 10, 3360 (2010).
  • Yang et al. (2010) M. Yang, A. Nurbawono, C. Zhang, and Y. P. Feng, Applied Physics Letters 96, 193115 (2010).
  • Balog et al. (2010) R. Balog, B. Jø rgensen, L. Nilsson, M. Andersen, E. Rienks, M. Bianchi, M. Fanetti, E. Laegsgaard, A. Baraldi, S. Lizzit, et al., Nature materials 9, 315 (2010).
  • Lin et al. (2015) C. Lin, Y. Feng, Y. Xiao, M. Dürr, X. Huang, X. Xu, R. Zhao, E. Wang, X.-Z. Li, and Z. Hu, Nano Lett. 15, 903 (2015).
  • Zhou et al. (2009) J. Zhou, Q. Wang, Q. Sun, X. S. Chen, Y. Kawazoe, and P. Jena, Nano letters 9, 3867 (2009).
  • Soriano et al. (2010) D. Soriano, F. Muñoz-Rojas, J. Fernández-Rossier, and J. J. Palacios, Phys. Rev. B 81, 165409 (2010).
  • Leconte et al. (2011) N. Leconte, D. Soriano, S. Roche, P. Ordejon, J.-C. Charlier, and J. J. Palacios, ACS Nano 5, 3987 (2011).
  • McCreary et al. (2012) K. M. McCreary, A. G. Swartz, W. Han, J. Fabian, and R. K. Kawakami, Physical Review Letters 109 186604(2012).
  • Soriano et al. (2011) D. Soriano, N. Leconte, P. Ordejón, J.-C. Charlier, J.-J. Palacios, and S. Roche, Phys. Rev. Lett. 107, 016602 (2011).
  • Nair et al. (2012) R. R. Nair, M. Sepioni, I.-L. Tsai, O. Lehtinen, J. Keinonen, A. V. Krasheninnikov, T. Thomson, A. K. Geim, and I. V. Grigorieva, Nature Physics 8, 199 (2012).
  • Saito et al. (1998) R. Saito, G. Dresselhaus, and S. Dresselhaus, Physical Properties of Carbon Nanotubes (Imperial College Press, 1998).
  • Yazyev and Helm (2007) O. V. Yazyev and L. Helm, Phys. Rev. B 75, 125408 (2007).
  • Boukhvalov et al. (2008) D. W. Boukhvalov, M. I. Katsnelson, and A. I. Lichtenstein, Phys. Rev. B 77, 035427 (2008).
  • Casolo et al. (2009) S. Casolo, O. M. Lø vvik, R. Martinazzo, and G. F. Tantardini, The Journal of chemical physics 130, 054704 (2009).
  • Yndurain (2014) F. Yndurain, Physical Review B 90, 245420 (2014).
  • Güttler et al. (2004) A. Güttler, T. Zecho, and J. Küppers, Surface Science 570, 218 (2004).
  • Palacios et al. (2008) J. J. Palacios, J. Fernández-Rossier, and L. Brey, Phys. Rev. B 77, 195428 (2008).
  • Jeloaica and Sidis (1999) L. Jeloaica and V. Sidis, Chemical Physics Letters 300, 157 (1999).
  • Hornekær et al. (2006) L. Hornekær, E. Rauls, W. Xu, Z. Sljivancanin, R. Otero, I. Stensgaard, E. Lægsgaard, B. Hammer, and F. Besenbacher, Phys. Rev. Lett. 97, 186102 (2006).
  • Chen et al. (2007) L. Chen, a. C. Cooper, G. P. Pez, and H. Cheng, Journal of Physical Chemistry C 111, 18995 (2007).
  • Denis and Iribarne (2009) P. a. Denis and F. Iribarne, Journal of Molecular Structure: THEOCHEM 907, 93 (2009).
  • Kerwin and Jackson (2008) J. Kerwin and B. Jackson, The Journal of chemical physics 128, 084702 (2008).
  • Ivanovskaya et al. (2010) V. V. Ivanovskaya, a. Zobelli, D. Teillet-Billy, N. Rougeau, V. Sidis, and P. R. Briddon, The European Physical Journal B 76, 481 (2010).
  • Sha and Jackson (2002) X. Sha and B. Jackson, Surface Science 496, 318 (2002).
  • Herrero and Ramirez (2009) C. P. Herrero and R. Ramirez, Physical Review B 79, 115429 (2009).
  • Huang et al. (2011a) L. F. Huang, M. Y. Ni, and Z. Zeng, J Phys Condens Matter 23, 435007 (2011).
  • Huang et al. (2011b) L. F. Huang, M. Y. Ni, G. R. Zhang, W. H. Zhou, Y. G. Li, X. H. Zheng, and Z. Zeng, J Chem Phys 135, 064705 (2011).
  • Borodin et al. (2011) V. A. Borodin, T. T. Vehviläinen, M. G. Ganchenkova, and R. M. Nieminen, Physical Review B 84, 075486 (2011).
  • Ferro et al. (2003) Y. Ferro, F. Marinelli, and A. Allouche, Chemical Physics Letters 368, 609 (2003).
  • Sljivancanin et al. (2009) Z. Sljivancanin, E. Rauls, L. Hornekaer, W. Xu, F. Besenbacher, and B. Hammer, The Journal of chemical physics 131, 084706 (2009).
  • Roman et al. (2007) T. Roman, W. A. Diño, H. Nakanishi, H. Kasai, T. Sugimoto, and K. Tange, Carbon 45, 218 (2007).
  • Andree et al. (2006) A. Andree, M. L. Lay, T. Zecho, and J. Küpper, Chemical Physics Letters 425, 99 (2006).
  • Lieb (1989) E. H. Lieb, Phys. Rev. Lett. 62, 1201 (1989).
  • Moaied et al. (2014) M. Moaied, J. V. Alvarez, and J. J. Palacios, Physical Review B 90, 115441(2014).
  • Hohenberg and Kohn (1964) P. Hohenberg and W. Kohn, Phys. Rev. 136, 864 (1964).
  • Kohn and Sham (1965) W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
  • Dion et al. (2004) M. Dion, H. Rydberg, E. Schröder, D. C. Langreth, and B. I. Lundqvist, Phys. Rev. Lett. 92, 246401 (2004).
  • Román-Pérez and Soler (2009) G. Román-Pérez and J. M. Soler, Phys. Rev. Lett. 103, 096102 (2009).
  • Kong et al. (2009) L. Kong, G. Román-Pérez, J. M. Soler, and D. C. Langreth, Phys. Rev. Lett. 103, 096103 (2009).
  • Ordejon et al. (1996) P. Ordejon, E. Artacho, and J. M. Soler, Phys. Rev. B 53, 10441 (1996).
  • Soler et al. (2002) J. M. Soler, E. Artacho, J. D. Gale, A. García, J. Junquera, P. Ordejon, and D. Sánchez-Portal, Journal of Physics: Condensed Matter 14, 2745 (2002).
  • Troullier and Martins (1991) N. Troullier and J. L. Martins, Phys. Rev. B 43, 1993 (1991).
  • Junquera et al. (2001) J. Junquera, Ó. Paz, D. Sánchez-Portal, and E. Artacho, Phys. Rev. B 64, 235111 (2001).
  • Zhang and Yang (1998) Y. Zhang and W. Yang, Physical Review Letters 80, 890 (1998).
  • Voter (2007) A. Voter, in Radiation Effects in Solids, edited by K. Sickafus, E. Kotomin, and B. Uberuaga (Springer Netherlands, 2007), vol. 235, chap. NATO Science Series, pp. 1–23.
  • Kalos and Whitlock (1986) M. H. Kalos and P. A. Whitlock, Monte Carlo Methods. Vol. 1: Basics (Wiley-Interscience, New York, NY, USA, 1986).
  • Bachellerie et al. (2007) D. Bachellerie, M. Sizun, D. Teillet-Billy, N. Rougeau, and V. Sidis, Chemical Physics Letters 448, 223 (2007).
  • Rougeau et al. (2006) N. Rougeau, D. Teillet-Billy, and V. Sidis, Chemical Physics Letters 431, 135 (2006).
  • Šljivančanin et al. (2012) Ž. Šljivančanin, R. Balog, and L. Hornekær, Chemical Physics Letters 541, 70 (2012).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description