Sticking of atomic hydrogen on graphene

Physisorption of hydrogen atoms on graphene and graphitic surfaces in general has long been considered a possible route for hydrogen storage. However, the depth of the physisorption well on clean surfaces is so small that no reliable storage device could be devised for operating at useful temperatures when only physisorbed species are considered (Tozzini and Pellegrini 2013). The problem remains of some interest for the chemistry of the interstellar clouds where formation of H₂ molecules on the interstellar grains may occur on low temperature surfaces where physisorbed species are stable (e.g. $T_{\rm s}=10$–20 K in the so-called diffuse clouds). At temperatures higher than $T_{\rm s}\sim 50$ K the rate of desorption is so large that refreshment of the surface is completed in a time-scale that is too short for astronomical standards.

The interaction potential relevant for this regime was characterized long ago by Ghio et al (1980) who used low energy H atom beams (50–65 meV) and first observed diffraction in the flux of scattered atoms off a graphite sample. They further measured a reduction of intensity in the specular reflection peaks, which was ascribed to the so-called selective adsorption, a kind of Feshbach resonance where the energy is temporarily channeled into closed diffraction channels. This confirmed the existence of an attractive interaction and suggested the presence of a reasonably deep adsorption well. From the position of the resonances Ghio et al estimated the well depth to be 43.3  ±  0.5 meV for graphite, and extrapolated this value to a single graphene layer to obtain 39.2  ±  0.5 meV (Ghio et al 1980).

Theoretical studies soon agreed in the position of the minimum being at the hollow site (i.e. at the center of a benzene ring) but the binding energies extracted from density functional theory (DFT) calculations were not sufficiently accurate because of well-known problems of DFT in handling dispersion forces. As a consequence, values in the range 0–100 meV were reported, despite the attempts of empirically correcting for van der Waals (vdW) interactions (see e.g. Ma et al (2011) and references therein). Good agreement between theory and experiment was achieved with the help of traditional quantum chemistry methods—second order, Møller–Plesset perturbation—on a cluster model (coronene, C₂₄H₁₂), using a rather large atomic basis-set and carefully correcting for the basis-set superposition error (BSSE) (Bonfanti et al 2007). This is essential to avoid sizable (unphysical) overbinding when weak interactions are of concern and geometry-dependent basis-sets are used. The H physisorption minimum was found at 2.93 Å above the surface plane with a depth of 39.7 meV at the hollow site, which decreases by only a few meV away from that site, thereby indicating a very small surface corrugation.

The findings of Bonfanti et al (2007) were confirmed by diffusion Monte Carlo studies on the same cluster model, although the same method applied to periodic graphene gave doubtful results (Ma et al 2011). Accurate parametrization of the empirical vdW correction was shown to give the required accuracy (Ferullo et al 2010) and, more recently, similar findings were obtained with the help of non-empirical, vdW-inclusive functionals (Bonfanti and Martinazzo 2018). These recent results pave the way for simulating H atom physisorption, diffusion and reactions with first principles means in a periodic setting, while accounting for the elusive vdW interactions in a reliable way.

Compared to several other surfaces, investigation of chemisorption of H atoms on graphitic substrates has started only recently (say the last 15 years), largely triggered by the raise of graphene. For some time, H atoms were not believed to be able to bind to the substrate, and early attempts to model the adsorption process without relaxing the surface failed in finding a chemisorption minimum. If fact, H atoms bind to the lattice only if (substantial) surface reconstruction is allowed, as first showed by Jeloiaca and Sidis on a cluster model (Jeloaica and Sidis 1999) and then by Sha and Jackson on a periodic setting (Sha and Jackson 2002). The authors of these studies found a binding energy $E_{\rm b}\sim 0.45$–0.65 eV, that was later refined with more converged calculations employing larger supercells. The value of the binding energy shows a sizable variability in the literature that can be ascribed to differences in the adopted DFT functionals and, more importantly, to the computational setup and the optimization strategy employed in the calculations. Accurate plane-waves DFT results unambiguously converge towards the value of $\sim$0.85 eV, (0.85 eV in a 4  ×  4 supercell (Hornekær et al 2006a), 0.84 eV in a 5  ×  5 supercell (Casolo et al 2009a) and 0.87 eV in a 8  ×  8 supercell (Lehtinen et al 2004)) while atomic-orbital based DFT results suggest a larger value (Ivanovskaya et al 2010) of 0.97 eV, likely because of the BSSE and the optimization strategy employed⁵. As for the length of the bond formed upon adsorption this is found typical of a CH in a hydrocarbon (∼1 Å), while the precise geometry of the deformed graphene structure remains sensitive to computational parameters (see table I in Casolo et al (2009a)).

The binding energy quoted above is of course much larger than the physisorbed one but remains rather small when compared to the interaction energy of H atoms on typical transition metal surfaces (∼2–2.5 eV) and to the CH bond energies in hydrocarbons (∼4 eV). The main reason for this discrepancy is that in forming a CH σ bond on graphene a considerable fraction of energy goes into the lattice, where it is stored indefinitely as deformation energy of the reconstructed (puckered) surface until the H atom is abstracted or desorbed, as we shall later when considering recombination of the adatoms to form H₂.

Such a surface reconstruction consists in a 0.3–0.4 Å out of plane displacement of the binding carbon atom, and occurs as a consequence of a sp²–sp³ re-hybridization of the valence C orbitals that is needed to form the CH bond. The re-hybridization induces a change in the geometry of the substrate site, from a planar (sp²) to a tetrahedral (sp³) form, thereby 'puckering' the surface plane, a distortion which actually extends for tens of Angstroms away from the adsorption site (see figure 1). This is also evident in the energy landscape of figure 2 which shows the energy as a function of the heights of both the H and the binding C atom, for a collinear configuration of the two atoms (the remaining degrees of freedom being left to relax): the deep chemisorbed well is found at a non-zero height of the carbon atom above the surface, and can thus be reached only if the C atom pulls out from the graphene plane.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Interestingly, the lattice distortion accompanying hydrogen chemisorption impacts also on the electronic properties of graphene. It greatly enhances the otherwise negligible spin–orbit coupling (Castro Neto et al 2009, Zhou et al 2010), and makes possible the realization of the spin Hall effect (Balakrishnan et al 2013).

Orbital re-hybridization is also responsible for an adsorption barrier  ∼0.2 eV high, as it is made evident in figure 2 where this barrier is seen to separate the chemisorbed region from the physisorbed one. Such barrier prevents direct H atom adsorption unless high-energy beams are employed, and allows physisorbed atoms to desorb rather than converting into stable chemisorbed species.

Neumann et al (1992) were among the first to realize the need of hyperthermal beams ($T_{\rm g} \sim 2300$ K) for depositing H atoms on natural graphite⁶. With the help of x-ray photoelectron spectroscopy (XPS) they unambiguously proved adsorption of H atoms on graphite with formation of a strong C-H bond. Further detailed information on hydrogen chemisorption were later provided by Zecho et al (2002a) when exposing highly oriented pyrolytic graphite (HOPG) to a high-flux hyperthermal beam ($T_{\rm g} \sim 2000$ K) of H(D) atoms, and using thermal desorption (TD) and high-resolution electron energy loss spectroscopy to investigate the species deposited on the surface at low temperature ($T_{\rm s}=150$ K). Puzzling at that time was the double peak structure in their TD spectra (Zecho et al 2002a) and the absence of any peak around $T_{\rm s}\sim 300$ K where desorption was expected according to the binding energy  ∼0.65 eV computed at that time. It took some years before dimer formation and recombination (Hornekær et al 2006a, 2006b) got uncovered and recognized to play a major role in the adsorption process (section 3).

From a theoretical perspective the barrier has been shown to arise from an avoided crossing between two diabatic electronic states, one purely repulsive in which the substrate electrons keep their ground-state coupling, and one strongly attractive describing the coupling between the H atom and a spin-excited electronic state of the substrate with two unpaired electrons (Casolo et al 2009a). The height of this barrier is yet unknown with precision, though results from standard DFT-GGA calculations cluster around 0.2 eV (Jeloaica and Sidis 1999, Sha and Jackson 2002, Ferro et al 2002, Hornekær et al 2006b, Casolo et al 2009a), a value that has been long considered a reasonable estimate. However, recent works employing vdW-inclusive DFT calculations claimed that the barrier height should be much smaller, if not vanishing (Moaied et al 2014, 2015, Brihuega and Yndurain 2018). Care is needed, though, in employing DFT for vdW interactions, especially in conjunction with atomic-like orbitals, because the combination of a (typically) overbinding functional with the ubiquitous BSSE may lead to errors of 0.2–0.3 eV that evidently wash out any adsorption barrier (Bonfanti and Martinazzo 2018). In addition, we mention that local and semilocal DFT functionals typically overestimate (underestimate) the binding (barrier) energies, particularly when the formation of true (i.e. directional) chemical bonds is involved, as it is the case for hydrogen on graphene. Coupled-cluster calculations on the coronene cluster model indeed found the barrier to be 370 meV, in reasonable agreement with previous multi-configuration quasi-degenerate perturbation theory results (Bonfanti et al 2011b). This estimate compares reasonably well only with the results of the accurate meta-hybrid functional M062X applied to the same system (328 meV), and it is considerable larger than the value obtained with the popular hybrid-functional B3LYP (262 meV) (Jensen et al 2018). Increasing the cluster size to circumcoronene (C₅₄H₁₈) decreases the barrier height to 304 (240) meV at the M062X (B3LYP) level, but it seems unlikely that its value attains  ∼0.2 eV when extrapolated to the infinite size limit. Interestingly, the meta-semilocal functional M06L, that has been finding increasing interest in the condensed matter community, provides even larger estimates for the barrier, 407 meV for coronene and 390 meV for circumcoronene (Jensen et al 2018).

From the experimental point of view, the height of the adsorption barrier is yet unknown, but its existence is undeniable. Without an energy barrier for adsorbing H atoms, cold atomic beams would be effective in depositing H atoms on a room temperature surface, something that has been ruled out with dedicated experiments using an inductively coupled plasma source delivering 0.025 eV H atoms (Aréou et al 2011). In addition, a barrierless adsorption would make hydrogen deposition a completely random process, with an estimated abundance of dimers (and larger clusters) much smaller than that observed experimentally. Rather, the existence of the barrier and its sensitivity to changes in the electronic structure of the substrate make the formation of certain adspecies configurations more likely than others, as we shall see in detail below, in section 3.

In closing this section we notice that chemisorption could be barrierless if C atoms were 'prepared' in a sp³ configuration, as it is partially the case in graphitic substrates that are keen to bend, e.g. crumpled single layer graphene, or that are already bent, e.g. carbon nanotubes. We shall discuss these issues below in section 5.

As mentioned above, a crucial step to investigate sticking of H atoms on graphene at energies comparable to or smaller than the barrier height, is to devise a tractable yet reliable model Hamiltonian describing dissipation into the lattice. This is possible upon exploiting the effective equivalence between the (generalized) Langevin dynamics and the Hamiltonian dynamics of a properly designed 'independent oscillator' (IO) model (Weiss 2008). For H sticking on graphene, the key assumption (that can be checked a posteriori, see for instance figure 3, left panel) is that the dynamics of the binding C atom can be described by a generalized Langevin equation (GLE),

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle m_{\rm C}\ddot{z}_{\rm C}(t)+m_{\rm C}\int_{-\infty}^{+\infty}\gamma(t-\tau)\dot{z}_{\rm C}(\tau){\rm d}\tau+\frac{\partial V(z_{\rm C}(t),{\bf x}_{\rm H}(t))}{\partial z_{\rm C}}=\xi(t)\label{eq:GLE} \nonumber \end{align} \tag{ 1 }$$

while the H atom is left to interact with C via an accurate (first principles) interaction potential $V$. In the above equation z_C is the height of the binding C (the degree of freedom most relevant for sticking), m_C its mass and ${\bf x}_{\rm H}$ is the position vector of the H atom. Furthermore, $\gamma(t)$ is understood to be a (causal) memory kernel and $\xi(t)$ is a zero-mean Gaussian stochastic process related to $\gamma(t)$ by a fluctuation-dissipation theorem of the second kind, $\renewcommand{\bra}[1]{\langle#1|} \renewcommand{\braket}[1]{\langle#1\rangle} \braket{\xi(t)\xi(0)}=\gamma(|t|)k_{B}T/m$. As consequence, the GLE above is entirely determined by the memory friction $\gamma(t)$ or, equivalently, by the so-called spectral density of the environmental coupling $J(\omega)$, that is defined to be $J(\omega)=m_{\rm C}\omega \Re \tilde{\gamma}(\omega)$, where $\tilde{\gamma}(\omega)$ is the Fourier transform of $\gamma(t)$ (Weiss 2008). Such spectral density turns out to be a key quantity for the non-Markovian Brownian dynamics of the GLE (Martinazzo et al 2011a, 2011b); in the present problem it can be understood as a phonon density of states on the binding carbon weighted by the coupling with the rest of the lattice.

Zoom In Zoom Out Reset image size

Equation (1) (and the accompanying equation for ${\bf x}_{\rm H}$) is clearly of classical nature. However, the above mentioned equivalence with an IO Hamiltonian dynamics makes it a valid starting point for designing a quantum Hamiltonian that governs the (dissipative) quantum dynamics of interest. For the present problem such Hamiltonian takes the form

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle H =\frac{{\bf p}_{\rm H}^{2}}{2m_{\rm H}}+\frac{p_{\rm C}^{2}}{2m_{\rm C}}+V({\bf x}_{\rm H},z_{\rm C})+\sum_{k=1}^{\rm F}\left[\frac{p_{k}^{2}}{2}+\frac{\omega_{k}^{2}}{2}\left(q_{k}-\frac{c_{k}}{\omega_{k}^{2}}(z_{\rm C}-z_{\rm C}^{\rm eq})\right)^{2}\right] \label{eq:IO final} \nonumber \end{align} \tag{ 2 }$$

and describes a hydrogen atom interacting with the binding carbon that, in turn, interacts with a set of oscillators mimicking the rest of lattice. Crucial for this description is the size of bath⁸ and the choice of the harmonic oscillator frequencies ($\omega_k$) and coupling coefficients (c_k). The latter are required to 'sample' the spectral density of the GLE, and are typically chosen according to a uniform frequency spacing, even though optimal sampling scheme can be easily devised (Martinazzo et al 2011b).

In turn, a successful modeling requires that the spectral density is derived from the atomistic description, and this can be accomplished with the help of classical, molecular dynamics simulations. In our work (Bonfanti et al 2015a) the starting point was the equilibrium autocorrelation function for the hydrogen displacement (or velocity) during its oscillatory motion in the chemisorbed state, since such information is in principle experimentally available through several vibrational spectroscopies. Furthermore, such choice does not require any a priori knowledge of the system potential (Bonfanti et al 2015a, 2015c, Bonfanti and Martinazzo 2016b, Gottwald et al 2016), hence suits well to an ab initio molecular dynamics (AIMD) approach. The resulting model, though derived from an equilibrium situation, proved to be robust enough to accurately describe the non-equilibrium sticking dynamics up to collision energies well above the height of the barrier (Bonfanti et al 2015b). This is shown in figure 3 (left panel) where the classical results obtained with the atomistic potential are compared with those obtained with the IO model.

The IO model of equation (2) was first used to study the vibrational relaxation of a H atom chemisorbed on graphene (Bonfanti et al 2015a). It was found that the surface mode describing block oscillations of the CH unit above the surface relaxes quickly, in few tens of fs, because its frequency (∼470 cm⁻¹) lies well below the Debye cutoff of the substrate. On the contrary, the H stretching mode (∼2550 cm⁻¹ at the PBE level, 2650 cm⁻¹ in the experiments, (Zecho et al 2002a)) lies above the Debye cutoff and thus its relaxation was found to be incomplete in the 3 ps-wide time window allowed by the chosen bath size (Bonfanti et al 2015a). Its lifetime was estimated to be  ∼5 ps, a value that is incidentally very close to the value of 5.2 ps obtained by Sakong and Kratzer applying Fermi's golden rule in a first principles setting (Sakong and Kratzer 2010), i.e. a more appropriate and accurate method for this weak-coupling regime.

The results of the quantum dynamical calculations (Bonfanti et al 2015b) are reported in figure 3, which shows the sticking coefficient for both the classical and the quantum case, as obtained in the limiting case in which the H atom impinges on top of the binding C (only for a $T_{\rm s}= 0$ K surface in the quantum case). Similar results have been recently obtained upon lifting the collinear assumption (Bonfanti et al 2018), apart for the size of the sticking probability that is approximately halved compared to the one in figure 3. In either cases it is found, as expected, that barrier-crossing plays the primary role at low collision energies, thereby allowing trapping into the chemisorption well for any projectile that is able to overcome the barrier. At high collision energies, however, energy transfer to the surface becomes the limiting factor, and fast H atoms hardly dissipate their excess energy and stick on the surface. As a consequence, the sticking coefficient attains a maximum at an energy which is about one and half larger than the barrier height.

Interestingly, it is further found that tunneling plays an important role only at energies much smaller than the barrier height. Rather, it is the zero-point motion of the binding carbon atom (i.e. the quantum fluctuations of the lattice) that plays a major role in determining sticking of the projectile H atoms, and makes a classical description of the process inadequate unless care is used to handle the zero-point energy issue. This is reasonable since, as figure 2 clearly shows, the motion of the C atom is strongly coupled to the 'sticking coordinate', and any tiny increase of its height above the surface determines a large decrease of the effective barrier to stick. This finding translates also in a marked dependence of the sticking probability on the surface temperature, but only in a classical setting: in the true, quantum setting the large Debye temperature of the surface leaves the lattice in the quantum regime for a wide temperature range.

As a matter of fact it was found that a simple impulsive model describing the collision of a classical projectile with a quantum surface reproduced the quantum results remarkably well for all but the lowest energies, thereby capturing the essential physics of such activated sticking dynamics, see right panel of figure 3 (Bonfanti et al 2015a, Bonfanti and Martinazzo 2016a). In this model, the energy available to the projectile for overcoming the barrier in the relative CH coordinate is the kinetic energy of H relative to C, and this is determined by both the collision energy E_coll and the thermal agitation of the binding surface atom. For the projectile atom to cross the barrier (assuming that it travels leftward towards C) the carbon atom speed v_C needs to exceed a threshold value $v_{\rm th}=-|v_{\rm H}|+v_{b}$, where $v_{\rm H}^{2}=2E_{\rm coll}/m_{\rm H}$ and $v_{b}^{2}=(1+\chi)2E_{\rm b}/m_{\rm H}$, $\chi=m_{\rm H}/m_{\rm C}$ being the projectile-target mass ratio and E_b the barrier height.

Once the atom has crossed the barrier, it is accelerated by the attractive interaction with the surface atom and energy transfer takes place to an extent that is determined by v_C, hence by the surface temperature T_s. The precise amount of energy that is transferred to the surface can be computed in the impulsive limit, which is appropriate at the high collision energies where energy dissipation dominates. Under these conditions, trapping occurs only when v_C lies in a specific interval

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle I_{i}(E_{\rm coll})=[v_{-},v_{+}],\ \ v_{\pm}(E_{\rm coll})=-\frac{1-\chi}{2}|\tilde{v}_{\rm H}|\pm\frac{1+\chi}{2}v_{0} \nonumber \end{align} \tag{ 3 }$$

where $\tilde{v}_{\rm H}^{2}=\frac{2(E_{\rm coll}+D)}{m_{\rm H}}$, $v_{0}^{2}=\frac{2(E_{\rm b}+D)}{m_{\rm H}}$ and D is an 'effective' depth for the interaction well. Hence, if this hard-collision limit kept down to low energies, the sticking probability would follow from integrating the distribution of the carbon atom velocities over the interval $\newcommand{\bi}{\boldsymbol}\Sigma(E_{\rm coll})=I_{i}(E_{\rm coll})\bigcap[v_{\rm th}(E_{\rm coll}), +\infty)$, i.e.

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle P_{s}(E_{\rm coll})=\int_{\Sigma(E_{\rm coll})}g(v){\rm d}v\label{eq:analytical P} \nonumber \end{align} \tag{ 4 }$$

where $g(v)$ is the velocity distribution appropriate to the C atom. Tunneling corrections are also possible at this stage by removing the threshold on v_C and introducing the appropriate tunneling probability in the integral, and this proved to improve the accuracy of the model down to the lowest energy considered (Bonfanti et al 2018).

Crucial for the success of the model is the choice of $g(v)$, since the classical (Maxwell–Bolztmann) velocity distribution proved to be inadequate to reproduce the quantum results. Rather, one needs the velocity distribution of the carbon atom quantum oscillator, which depends on both its fundamental frequency and on the coupling with the rest of the lattice. Since v_C is a linear combination of Gaussian variables (the normal modes of the substrate) the required function $g(v)$ is readily seen to take a Gaussian form similar to that of an independent oscillator (Bonfanti et al 2015b),

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle g(v)=\sqrt{\frac{m_{c}}{\pi\hbar\Omega_{T}}}{\rm e}^{-\frac{m_{c}v^{2}}{\hbar\Omega_{T}}}\label{v-distribution} \nonumber \end{align} \tag{ 5 }$$

but with a temperature-dependent effective frequency $\Omega_{T}$ that accounts for the coupling to the bath. Only in the high temperature limit (in fact, for $T_{\rm s}\gg T_{\rm D}/2$, where T_D is the Debye temperature, $T_{\rm D}\sim400$ K in graphene) $\Omega_T$ increases linearly with T_s and $g(v)$ in equation (5) takes the standard form of a Maxwell–Boltzmann distribution.

The computed 'initial' sticking co-efficient has been recently used in a simple kinetic model devised to understand the adsorption process of H atoms on graphene at low coverages (Bonfanti et al 2018). The model included only sticking, dimer formation and Eley–Rideal reactions, and was found to describe rather well the hydrogen uptake curves measured long ago by Beitel on graphite (Beitel 1969) and more recently by Haberer and coworkers on graphene on gold (Haberer et al 2011a, 2011b). It also helped to rationalize the isotopic effect observed by Paris et al (2013) but, unexpectedly, was found somewhat at odds with the findings on HOPG (Zecho et al 2002a, Andree et al 2006). The reason of this discrepancy might lie in the fact that HOPG is of high-quality form only from a crystallographic point of view, and its defect concentration might have a pronounced effect on the macroscopic adsorption process.

For hydrogen recombination on graphene LH is relevant in the physisorbed regime only since, as shown above, chemisorbed atoms prefer to desorb rather than diffusing. Even in this case, however, LH is not really standard, since thermalization of the physisorbed species in stable adsorption sites (as required by the reaction mechanism) is hampered by zero-point fluctuations. The reaction efficiency for two atoms trapped on the surface and approaching each other has been studied with wave packet techniques (Morisset et al 2004a, 2005), and found very high. The reaction should occur through the deflection of the two projectiles toward and away from the surface plane: the rebound of the first transfers energy to the nascent molecule giving rise to a strong rotational excitation. Then, molecular hydrogen either desorbs immediately or stays trapped on the surface in a metastable state, and it is vibrationally hot, and rotationally highly excited (Morisset et al 2004a, 2005).

The relevance of HA recombination for hydrogen on graphene is yet unclear, for it requires a weakly corrugated but rather strong atom–surface interaction that can trap hyperthermal atoms and prevent them to desorb during their excursion along the surface. HA is probably operative when H atoms trap in the physisorbed well of well-cleaned surfaces, but it is expected to be rather sensitive to the surface temperature because of the huge effect that the latter has on the hot-atom lifetime. In addition, formation of hot-atoms from gas-phase projectiles requires a source of corrugation (such as, e.g. that provided by an adsorbed species) that can channel energy from the beam direction (typically normal to the surface) to the surface plane. Nevertheless, there are evidences that HA reaction occurred, in conjunction with ER abstraction, in the experiments by Zecho et al on HOPG (Zecho et al 2002b), since the extracted 'effective' abstraction cross-sections (>16 $\mathring{\rm A}^2$ at $T_{\rm s}=150$ K) were found to decrease with increasing surface temperature and to level off at high temperature (∼3.5 $\mathring{\rm A}^2$) (Zecho 2007). These results suggest that hot-atoms might play an important role when the temperature is sufficiently low (and the surface sufficiently clean) that such fast-moving atoms stay on the surface long enough to encounter a reaction partner.

One comment is in order in this context since the same abstraction experiments (Zecho et al 2002a) were found to agree rather well with the results of the first realistic quantum calculations of the ER cross-section (Zecho et al 2002b). The agreement should be considered fortuitous since, on one hand, the theoretical results were not corrected for the 1/4 spin-statistical factor appropriate for this reaction (Casolo et al 2009b, 2013) (see also below) while, on the other hand, as argued above, the experimental ones likely accounted for HA reactions too. Comparison between the saturation value of the cross-section (∼3.5 $\mathring{\rm A}^2$) and 1/4 of the theoretical estimate would give an agreement similar to that claimed in the original work (Zecho et al 2002b), which can be improved by accounting for the competing reaction channels (Casolo et al 2013).

The ER reaction, usually of secondary importance for different surfaces, becomes very important here for chemisorbed species, in the wide temperature range where the latter are stable on the surface and physisorbed species are not present (say for T_s in the range 50–300 K at low coverage, and T_s in the range 50–500 K at higher coverage, see section 3). As will be discussed in section 3.1.2, it generally competes with sticking (dimer formation) when the substrate is exposed to hot H atoms, and the 'branching ratio' depends markedly on the energy of the incoming H atoms. However, when cold atoms are used on a H pre-covered sample, abstraction dominates over dimer formation, as shown experimentally (Aréou et al 2011) and confirmed by theory (Casolo et al 2013).

This reaction mechanism has been extensively studied over the years, both theoretically and experimentally, mainly because of its importance for the interstellar chemistry (Tielens 2013). Several specific aspects of the dynamics have been addressed (the size of the cross sections, the internal excitation of the product molecules, the role of the collision energy and of the vibrational excitation of the adsorbate, the effect of isotopic substitutions, etc) upon resorting to various approximations, either in the dynamics or in the model (Jackson and Lemoine 2001, Meijer et al 2001, Farebrother et al 2002, Zecho et al 2002b, Sha et al 2002, Morisset et al 2003, 2004b, Martinazzo and Tantardini 2006a, 2006b, Casolo et al 2009b, Bachellerie et al 2009b, Sizun et al 2010, Pasquini et al 2016). Quantum dynamics was mostly restricted to the rigid, flat approximation (Sha et al 2002, Zecho et al 2002b, Martinazzo and Tantardini 2005, 2006a, Casolo et al 2009b), and the dynamics was addressed in two opposite limits, namely a reaction much faster than the surface atom motion ('diabatic' limit) or such slow that C fully relaxes during the projectile motion ('adiabatic' limit). These quantum studies all show a reaction cross-section that increase steadily with energy till the competing collision induced desorption (CID) channel opens up. In this regime, total reaction cross sections showed singular quantum oscillations that witnessed an unusual reaction mechanism leading to selective population of the low-lying H₂ vibrational levels (Martinazzo and Tantardini 2005, 2006a). Extension to the cold collision energy regime $E_{\rm coll}\sim10^{-4}$–10⁻² eV (Casolo et al 2009b), on the other hand, showed that reaction could also be hampered by quantum reflection despite the absence of any barrier (Bonfanti et al 2011a). This quantum effect may indeed arise in the presence of deep and narrow potential wells whenever the projectile De Broglie's wavelength approaches the size of the well to be crossed.

Overall, general consensus has been reached on the large size of the cross-section and on the flow of a large fraction of the reaction exothermicity into product vibrational excitation. The strong H–H interaction dominates the dynamics, and product molecules can form and leave the surface as soon as the two hydrogen atoms get closer than about twice the equilibrium internuclear distance in H₂ (∼0.7 Å). In addition, chemisorbed H atoms are relatively far from the surface and this allows steering of the projectile, with relatively large cross-section (∼10 $\mathring{\rm A}^2$ when not yet corrected for the spin statistical factor). This is in sharp contrast with metal surfaces where stronger H–substrate interactions (and shorter equilibrium geometries) 'mask' the target atom to the hydrogen projectile and leads to much smaller cross-section (<1 $\mathring{\rm A}^2$ on metals).

For physisorbed targets the reaction was found to be very efficient down to 1 K (Martinazzo and Tantardini 2006b, Casolo et al 2009b). However, in this case the CID channel opens up at quite small energies (∼40 meV) and becomes soon very efficient, thereby causing a rapid drop of the reaction efficiency. A large cross-section for trapping H atom projectiles was predicted that might trigger HA reactions (Martinazzo and Tantardini 2006b, Casolo et al 2009b).

As mentioned above, a striking feature of the H₂ formation reactions on graphene and graphitic surfaces is the vibrational excitation of the product molecules. This is interesting for several reasons. In the chemistry of the interstellar clouds, for instance, the vibrational excitation of the H₂ molecules triggers some endothermic reactions that would be otherwise prohibited by the severe conditions of the interstellar clouds, e.g. H₂  +  C⁺   →  CH⁺   +  H (Agúndez et al 2010, Bonfanti et al 2014). On the other hand, vibrationally hot molecules easily form negative ions by dissociative attachment, H₂  +  e⁻  →  H  +  H⁻. Thus, they provide a viable route to produce negative ion sources, to be used, for instance, in heating and current drive systems in experimental fusion reactors.

Theoretical calculations agree on the vibrational excitation of the product molecules, being them obtained from LH or ER, employing either chemisorbed or physisorbed species, although the precise level of excitation may differ from one study to the other because of details in the interaction potentials or in the adopted dynamical approach. Experiments employed cold HOPG surfaces ($T_{\rm s} = 10$–15 K) and cold beams ($T_{\rm g} \sim 300$ K), i.e. in a regime mostly relevant for LH recombination between physisorbed species, and used resonant enhanced multi photon ionization to probe the product energy (Latimer et al 2008). The rovibrational distributions of the nascent molecule were found to peak at the first few values of the rotational quantum number j and at much larger value of the vibrational quantum number $\nu=4, 5$. The observed low rotational excitation is at odds with theoretical predictions on LH/HA recombination of physisorbed species, since in the simulations H₂ molecules were found with a much larger rotational excitation (Morisset et al 2004a, 2005, Bachellerie et al 2009a). They agree much better with those computed with ab initio molecular dynamics in an unrestricted investigation of the ER reaction involving chemisorbed species (Casolo et al 2013), thereby suggesting that chemisorbed H atoms might play also a role at the extreme conditions of the experiments (Latimer et al 2008), see figure 4 and section 3.1.2.

Zoom In Zoom Out Reset image size

The dynamical role of the lattice, as well as its ability to absorb part of the reaction energy, has received little attention. Molecular dynamics (Bachellerie et al 2009a) and ab initio molecular dynamics (Casolo et al 2013, 2016) investigations did assess the role of the substrate but in a classical setting and with different aims, while, as mentioned above, quantum dynamical studies were mostly restricted to reduced dimensional models that completely neglected the dynamical role of the substrate or reduced it to that of the binding carbon only at the expense of other degrees of freedom.

Progress in this direction has been recently made by extending the model developed for the sticking dynamics (section 2.2) to account for the presence of an additional H atom (Pasquini et al 2018). Despite the limitation of a collinear approach, the mentioned study represents the first unrestricted investigation of the C atom dynamics (and its energy exchange with the rest of the lattice) during the reaction, in the correct quantum setting most appropriate for a light atom dynamics. Indeed, even though the reaction dynamics is reasonably well described by classical means at all but the lowest collision energies, lattice quantum fluctuations need to be correctly taken into account for a correct description of the reaction (see also section 2.2). One of the main findings of this study is that the C atom dynamics is essential for the reaction, while the rest of the lattice plays a secondary role only. The reason is that, even though unpuckering of the surface is fast—some tens of fs, in accordance with the short relaxation time found for the related CH surface mode, see section 2.2—it starts only after formation of H₂ is completed. As a major consequence a considerable substrate heating occurs: the energy stored in the puckered surface is left entirely on the lattice, ∼0.8 eV per reactive event. Adding the energy that is necessarily dissipated when chemisorbing the first H atom (the 'target' in the ER abstraction) one finds that about 1.6 eV per H atom pair is stored in graphene. This result is in sharp contrast with the situation in which two physisorbed H atoms recombine via LH kinetics, in which case only (twice) the H atom physisorption energy is left on the surface.

It is worth noticing in this context that, although theoretical modeling was invariably performed in the (electronically) adiabatic approximation, it is likely that the large amount of chemical energy that is quickly released into the lattice triggers non-adiabatic effects and creates e–h pairs. In principle, such 'chemically induced' e–h pairs could be detected as chemi-currents (Gergen 2001, Nienhaus 2002, Kandratsenka et al 2018).

From a chemical perspective midgap states arise naturally in the resonant valence bond (RVB) description of graphene once breaking a double bond and 'resonating'. Such a picture—though most often confined to a qualitative level—is a simple graphical translation of the RVB variational wavefunction which, for a singlet state, reads as

$$\begin{align*} \newcommand{\e}{{\rm e}} \renewcommand{\ket}[1]{|#1\rangle} \displaystyle \ket{\Psi}=\sum_{I}c_{I}\Pi_{ij}^{\{I\}}\left(a_{i\uparrow}^{\dagger}a_{j\downarrow}^{\dagger}+a_{j\uparrow}^{\dagger}a_{i\downarrow}^{\dagger}\right)\ket{0} \nonumber \end{align*}$$

where $a_{i\sigma}^{\dagger}$ creates an electron in site i with spin σ, the product runs over all pairs ($i\neq j$) that define linearly independent 'spin-couplings' and c_I's are variational parameters. In the above equation each pairing I represents a chemical structure and the superposition is commonly known as chemical resonance. Figure 6, for instance, shows two of the many graphene's chemical structures with $\leftrightarrow$ indicating the above superposition. The number of independent couplings is a group-theoretical object that increases quickly with increasing the number of electrons¹⁰ and this prevents the use of the above ansatz for all but the simplest systems. However, limiting the pairs to the nearest neighbors ones (the so-called Kekulé structures) leads to considerable saving without seriously affecting the accuracy (Wu et al 2002, 2003) (coronene, for instance, with its 24 π electrons has 208,012 singlet couplings but only 20 Kekulé structures) and further saving is possible by recognizing that the largest contribution to the chemical resonance is provided by the Clar's formulas (Solà 2013), i.e. those structures that maximize the number of sextets (for coronene there exist only two of such Clar's formulas).

Zoom In Zoom Out Reset image size

Now, considering H adsorption, it is not hard to realize that binding of a H atom breaks the aromatic network and leaves one unpaired electron on the lattice that is free to move by 'bond switching': spin re-coupling with a neighboring double bond creates an unpaired electron in one every two lattice sites, on the same majority set predicted by the counting rule mentioned above (figure 6, bottom). In fact, such itinerant electron is equivalently described by a singly occupied molecular orbital delocalized on such a set. As we shall see in the following, despite its simplicity, such picture is powerful in describing spin and/or charge density localization in π-conjugated molecules and in graphene, and proved to be rather useful for predicting chemical reactivity.

A quantitative look at the zero-energy modes in graphene is provided by the (noninteracting) Anderson model (Anderson 1961, Robinson et al 2008, Wehling et al 2010), which we discuss here in some detail because it is instrumental for the following section. In this model the ad-species is described by a single level at energy $\newcommand{\e}{{\rm e}} \epsilon_{{\rm ad}}$ (in our case the 1 s level of the H atom), which is let to hybridize with a carbon atom of the lattice, say of A type at the origin, by adding the term

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle H_{{\rm ad}}=\sum_{\sigma}\epsilon_{{\rm ad}}d_{\sigma}^{\dagger}d_{\sigma}+W\sum_{\sigma}\left(a_{0,\sigma}^{\dagger}d_{\sigma}+d_{\sigma}^{\dagger}a_{0,\sigma}\right) \nonumber \end{align*}$$

to the lattice Hamiltonian

$$\begin{align*} \newcommand{\e}{{\rm e}} \renewcommand{\braket}[1]{\langle#1\rangle} \renewcommand{\bra}[1]{\langle#1|} \displaystyle H_{{\rm latt}}=-t\sum_{\sigma}\sum_{\braket{i,j}}a_{i,\sigma}^{\dagger}a_{j,\sigma}+t'\sum_{\sigma}\sum_{\braket{\braket{i,j}}}a_{i,\sigma}^{\dagger}a_{j,\sigma}+... \nonumber \end{align*}$$

In these equations W is the hybridization energy, $d_{\sigma}^{\dagger}$ $(d_{\sigma})$ creates (destroys) an electron with spin σ in the adatom energy level, $a_{i, \sigma}^{\dagger}$ ($a_{i, \sigma}$) does the same for the lattice site i and t, $t'$ are the hopping energies for the nearest-neighbors (NN) and next-nearest-neighbors pairs. The problem is best handled with the help of the Green's operator¹¹ $G(\lambda)$ of the one-electron Hamiltonian, upon partitioning the one-electron space into a primary subspace (the lattice) and the remainder (Levine 1969). In fact, it is not hard to show that the projection of the exact Green's operator onto the lattice takes the simple form¹² $G_{{\rm ef\,\!f}}(\lambda)=(\lambda-H_{{\rm ef\,\!f}}){}^{-1}$ where $H_{{\rm ef\,\!f}}$ is an effective, energy-dependent Hamiltonian implicitly accounting for the dynamics in the adatom level, i.e. $H_{{\rm ef\,\!f}}=H_{{\rm latt}}+V(\lambda)$, where $V(\lambda)$ reads as

$$\begin{align*} \newcommand{\e}{{\rm e}} \renewcommand{\ket}[1]{|#1\rangle} \renewcommand{\bra}[1]{\langle#1|} \displaystyle V(\lambda)=W^{2}\sum_{\sigma}\frac{\ket{\chi_{0,\sigma}}\bra{\chi_{0,\sigma}}}{\lambda-\epsilon_{{\rm ad}}} \nonumber \end{align*}$$

$\renewcommand{\ket}[1]{|#1\rangle} \ket{\chi_{0, \sigma}}$ being a p_Z spin–orbital at site 0. It is thus seen that, as long as the electron dynamics in the lattice is of concern, hybridization with an impurity level introduces an energy dependent scattering potential $V$ that becomes strong for energies approaching the adatom level. A solution of this scattering problem is possible with the help of the T operator (Taylor 1969), here defined to be $T(\lambda)=V+VG_{{\rm ef\,\!f}}(\lambda)V$, since the corresponding Lippmann–Schwinger equation $T(\lambda)=V+VG^{0}(\lambda)T(\lambda)$ (where G⁰ is the Green's operator of the unperturbed lattice) is solved by

$$\begin{align} \newcommand{\e}{{\rm e}} \renewcommand{\ket}[1]{|#1\rangle} \renewcommand{\bra}[1]{\langle#1|} \displaystyle T(\lambda)=t(\lambda)\sum_{\sigma}\ket{\chi_{0,\sigma}}\bra{\chi_{0,\sigma}}\ \ \ \ t(\lambda)=\frac{W^{2}}{\lambda-\epsilon_{{\rm ad}}-W^{2}g_{00}^{0}(\lambda)}\label{eq:T-matrix} \nonumber \end{align} \tag{ 6 }$$

when the potential takes a separable form. In turn, knowledge of the T-matrix allows one to obtain the required lattice Green's operator as $G_{{\rm ef\,\!f}}(\lambda)=G^{0}(\lambda)+G^{0}(\lambda)T(\lambda)G^{0}(\lambda)$ and thus solve the scattering problem. In equation (6), $g_{00}^{0}(\lambda)$ is the on-site Green's function of the unperturbed lattice, $\renewcommand{\bra}[1]{\langle#1|} \renewcommand{\braket}[1]{\langle#1\rangle} g_{00}^{0}=\braket{{\bf 0}\sigma|G^{0}(\lambda)|{\bf 0}\sigma}$, which for graphene in the linear band approximation takes the form

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle g_{00}^{0}(\epsilon)\approx-\frac{\epsilon}{\epsilon_{c}^{2}}\ln\left(\left|\frac{\epsilon_{c}^{2}}{\epsilon^{2}}-1\right|\right)-{\rm i}\pi\rho^{0}(\epsilon) \nonumber \end{align*}$$

where $\newcommand{\e}{{\rm e}} \rho^{0}(\epsilon)$ is the density of states per C atom per spin channel, $\newcommand{\e}{{\rm e}} \rho^{0}(\epsilon)=\frac{|\epsilon|}{\epsilon_{c}^{2}}\Theta(\epsilon_{c}-|\epsilon|)$. Here, $\newcommand{\e}{{\rm e}} \epsilon_{c}$ is an energy cutoff that determines the bandwidth and that is inversely related to the carbon–carbon bond length d_CC, i.e. $\newcommand{\e}{{\rm e}} \epsilon_{c}\approx\hbar v_{\rm F}d_{\rm CC}^{-1}$ (matching $\newcommand{\e}{{\rm e}} \rho^{0}(\epsilon)$ to the known low-energy expression for the density of states in graphene would result in $\newcommand{\e}{{\rm e}} k_{c}=\epsilon_{c}/\hbar v_{\rm F}=2\sqrt{2\pi/3\sqrt{3}}d_{\rm CC}^{-1}\approx2d_{\rm CC}^{-1}$ or $\newcommand{\e}{{\rm e}} \epsilon_{c}=t\sqrt{\pi\sqrt{3}}\approx6$ eV). Notice that for $\newcommand{\e}{{\rm e}} \epsilon\rightarrow 0$, for any linear (and sublinear) density of states the real part of the on-site Green's function ($\newcommand{\e}{{\rm e}} \Re g_{00}^{0}(\epsilon)$) features a vertical cusp that makes the position of the resonance in the T-matrix (i.e. the lowest energy solution of the equation $\newcommand{\e}{{\rm e}} \epsilon-\epsilon_{{\rm ad}}-W^{2}\Re g_{00}^{0}(\epsilon)=0$) rather insensitive to the hybridization energy W and always very close to $\newcommand{\e}{{\rm e}} \epsilon=0$.

With the $G(\lambda)$ operator at hand, one readily obtains the density of states at the lattice sites, $\newcommand{\e}{{\rm e}} \rho_{i, \sigma}(\epsilon)=-\frac{1}{\pi}\Im G_{i\sigma, i\sigma}(\epsilon)$. The calculation is straightforward when the total density of states is of interest, since elementary algebra shows that it takes the form (per C atom, per spin channel)

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \rho(\epsilon)=\rho^{0}(\epsilon)+\frac{1}{N_{\rm C}}\frac{1}{\pi}\Im\left[\frac{\partial{\rm ln}t(\epsilon)}{\partial\epsilon}\right] \nonumber \end{align*}$$

where N_C is the number of lattice sites (this expression contains also the contribution of the adatom level, $\newcommand{\e}{{\rm e}} \rho_{d, \sigma}(\epsilon)=-\frac{1}{\pi}\Im G_{d\sigma, d\sigma}(\epsilon)$, as it turns out from the Green's operator projected in the secondary space, $\newcommand{\e}{{\rm e}} G_{d\sigma, d\sigma}(\lambda)=1/\left(\lambda-\epsilon_{{\rm ad}}-W^{2}g_{00}^{0}(\lambda)\right)\equiv t(\lambda)/W^{2}$). It is readily seen that $\newcommand{\e}{{\rm e}} \Delta n(\epsilon)=\rho(\epsilon)-\rho^{0}(\epsilon)<0$ when $\newcommand{\e}{{\rm e}} \epsilon\approx0$, except for a positive peak $\newcommand{\e}{{\rm e}} \sim\epsilon_{c}^{2}/2W^{2}|\epsilon_{{\rm ad}}|$ wide around the adatom energy level when $\newcommand{\e}{{\rm e}} W<\epsilon_{c}/\sqrt{2}$ which becomes a sharp peak between $\newcommand{\e}{{\rm e}} \epsilon=0$ and $\newcommand{\e}{{\rm e}} \epsilon={\rm sign}(\epsilon)\epsilon_{c}^{2}/2W^{2}$ when $\newcommand{\e}{{\rm e}} W>\epsilon_{c}/\sqrt{2}$. Wehling et al (2009) fitted ab initio results to tight binding models and found for hydrogen (and several other neutral species) $\newcommand{\e}{{\rm e}} \epsilon_{{\rm ad}}\approx-0.2$ eV and $W\approx2t=5.2$ eV, thereby placing the resonance at a energy of about  −0.06 eV, in agreement with the results of figure 5. In this context, it is worth noticing that both the adatom level and the binding C site contribute little to the density of states in this spectral window, since they form a bonding–antibonding pair of states which move far from the Fermi level. This becomes evident when $\newcommand{\e}{{\rm e}} W\gg\epsilon_{c}$ since in such case the binding C and the adatom level give each a contribution of $\newcommand{\e}{{\rm e}} \frac{1}{2}\left[\delta(\epsilon-W)+\delta(\epsilon+W)\right]$ which separates out from the band.

Of great interest are also the spatial properties of the H-induced resonance which, as seen in figure 5, semilocalizes around the adatom and decays away from it with a characteristic power law. This behavior is also evident when imaging H atoms with STM, since a bright protusion around the adatom appears at small bias and displays a characteristic threefold symmetry (Wang et al 2006). In the T-matrix approach to the Anderson problem these properties can be investigated by looking at the scattered component $\renewcommand{\ket}[1]{|#1\rangle} \newcommand{\e}{{\rm e}} \ket{\psi_{\epsilon}^{{\rm scatt}}}$ of the scattering state $\renewcommand{\ket}[1]{|#1\rangle} \newcommand{\e}{{\rm e}} \ket{\psi_{\epsilon}+}$,

$$\begin{align*} \newcommand{\e}{{\rm e}} \renewcommand{\ket}[1]{|#1\rangle} \displaystyle \ket{\psi_{\epsilon}+}=\ket{\psi_{\epsilon}^{0}}+\ket{\psi_{\epsilon}^{{\rm scatt}}}=\ket{\psi_{\epsilon}^{0}}+G^{0}(\epsilon)T(\epsilon)\ket{\psi_{\epsilon}^{0}} \nonumber \end{align*}$$

where $\renewcommand{\ket}[1]{|#1\rangle} \newcommand{\e}{{\rm e}} \ket{\psi_{\epsilon}^{0}}$ is an eigenstate of the unperturbed lattice with energy . It is clear that such scattering component requires computation of the off-site Green's function, being it proportional to

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \psi_{\epsilon}^{{\rm scatt}}({\bf r})\propto t(\epsilon)G^{0}({\bf r},{\bf 0}|\epsilon) \nonumber \end{align*}$$

where $\renewcommand{\bra}[1]{\langle#1|} \renewcommand{\braket}[1]{\langle#1\rangle} \newcommand{\e}{{\rm e}} G^{0}({\bf r}, {\bf 0}|\epsilon)=\braket{{\bf r}\sigma|G^{0}(\epsilon)|{\bf 0}\sigma}$ for a H adatom binding to the site at the origin (of A type in the following). This quantity has been considered by several authors for different reasons, ranging from investigations of STM imaging (Wang et al 2006) to Ruderman–Kittel–Kasuya–Yosida (RKKY) interactions in graphene (Sherafati and Satpathy 2011a, 2011b, Kogan 2011). It can be obtained numerically as Fourier transform of the (simpler) Green's function in k-space, namely from

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle G_{XA}^{0}({\bf r},{\bf 0}|\epsilon)=\frac{1}{\Omega_{\rm BZ}}\int_{\rm BZ}{\rm d}^{2}{\bf k}{\rm e}^{{\rm i}{\bf k}{\bf r}}G_{XA}^{0}({\bf k}|\epsilon) \nonumber \end{align*}$$

where $X=A, B$ depending on whether ${\bf r}$ is a lattice position in the A or B sublattice ($\mathcal{A}$ and $\mathcal{B}$, respectively, in the following), the integral runs over the Brillouin zone (BZ) and $\Omega_{\rm BZ}$ is its area. Here, $\renewcommand{\bra}[1]{\langle#1|} \renewcommand{\braket}[1]{\langle#1\rangle} \newcommand{\e}{{\rm e}} G_{XY}^{0}({\bf k}|\epsilon)=\lim_{\eta\rightarrow0^{+}}\braket{\phi_{{\bf k}X}|(\epsilon+{\rm i}\eta-H_{{\bf k}}){}^{-1}|\phi_{{\bf k}Y}}$, where $\renewcommand{\ket}[1]{|#1\rangle} \ket{\phi_{{\bf k}, X}}$ is a Bloch state built with $X=A, B$ basis functions and $H_{{{\rm \bf k}}}$ is the k  −  Hamiltonian (for the above two-band model of graphene, a 2  ×  2 matrix in the sublattice indexes). Fortunately, at the low energies of interest for the H-induced resonance, one can use the linear band approximation and expand the above integrand around the K, K' points, remove the energy cutoff $\newcommand{\e}{{\rm e}} \epsilon_{c}$ with negligible error (for not too small distances, $r>d_{\rm CC}$) and perform the integral analytically. Following this route we find

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle G_{AA}^{0}({\bf r},{\bf 0}|\epsilon)=-{\rm i}\frac{A_{c}}{2}\frac{|\epsilon|\cos({\bf K}{\bf r})}{\hbar^{2}v_{\rm F}^{2}}H_{0}^{\pm}\left(\frac{|\epsilon|r}{\hbar v_{\rm F}}\right)\ \ \ \ {\rm for}\ {\bf r}\in\mathcal{A} \nonumber \end{align*}$$

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle G_{BA}^{0}({\bf r},{\bf 0}|\epsilon)=+{\rm i}\frac{A_{c}}{2}\frac{\epsilon\sin({\bf K}{\bf r}+\theta)}{\hbar^{2}v_{\rm F}^{2}}H_{1}^{\pm}\left(\frac{|\epsilon|r}{\hbar v_{\rm F}}\right)\ \ \ \ {\rm for}\ {\bf r}\in\mathcal{B} \nonumber \end{align*}$$

where $H_{l}^{\pm}$ are Hankel functions of the first and second kind (for positive and negative energies, respectively) of order l, A_c is the area of the graphene unit cell, θ is the angle between ${\bf r}$ and ${\bf K}$ and, finally, ${\bf K}$ locates a K corner in the BZ. The latter can be given as $2\pi{\bf K}=\Omega_{\rm BZ}\boldsymbol{\delta}\wedge\hat{{\bf n}}$, where $\boldsymbol{\delta}$ is the position of the B site relative to the A one in the (arbitrarily) chosen unit cell, and $\hat{{\bf n}}$ is the surface normal.

One thus sees that the scattering resonance has the required threefold symmetry, with maxima in the three directions orthogonal to the ${\bf K}$'s equivalent corners of the BZ, i.e. along the armchair directions¹³. Importantly, in the intermediate-distance range $\newcommand{\e}{{\rm e}} d_{\rm CC}<r\ll\hbar v_{\rm F}/|\epsilon|$ (which can be rather wide at the energies of interest for an impurity-induced resonance) one can replace the Hankel's functions with their low-argument expansion obtaining

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle G_{AA}^{0}({\bf r},{\bf 0}|\epsilon)\approx\frac{A_{c}}{\pi}\frac{\epsilon r}{\hbar v_{\rm F}}\frac{\cos({\bf K}{\bf r})}{\hbar v_{\rm F}r}\left[\ln\left(\frac{|\epsilon|r}{2\hbar v_{\rm F}}\right)+\gamma\right]\ \ \ \ {\rm for}\ {\bf r}\in\mathcal{A} \nonumber \end{align*}$$

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle G_{BA}^{0}({\bf r},{\bf 0}|\epsilon)\approx\frac{A_{c}}{\pi}\frac{\sin({\bf K}{\bf r}+\theta)}{\hbar v_{\rm F}r}\ \ \ \ {\rm for}\ {\bf r}\in\mathcal{B} \nonumber \end{align*}$$

where $\gamma\approx0.577$ is the Eulero–Mascheroni constant. The striking feature of this expression is the 1/r decay of G⁰ (hence of the resonance wavefunction $\newcommand{\e}{{\rm e}} \psi_{\epsilon}^{{\rm scatt}}$) on the B sublattice and its smaller magnitude on the A sublattice that hosts the adatom (with $\newcommand{\e}{{\rm e}} x=|\epsilon|r/\hbar v_{\rm F}$ one has $|x\ln x|\leqslant e^{-1}$ for x  <  1), vanishing for $\newcommand{\e}{{\rm e}} \epsilon\rightarrow0$. Such algebraic decay of the wavefuncion was first predicted (Cheianov and Fal'ko 2006, Pereira et al 2006, Mariani et al 2007, Pereira et al 2008, Bena 2008) and experimentally observed (Ugeda et al 2010) for carbon atom vacancies. Detailed experiments on hydrogen atoms have been recently performed by Gonzalez-Herrero et al (2016), who observed and manipulated H atoms adsorbed on graphene with scanning tunneling microscopy (STM)/scanning tunneling spectroscopy techniques. Their measured differential conductance maps around the adatoms provide a clear representation of the spatial localization of the impurity-induced state, as is clearly seen in figure 7.

Zoom In Zoom Out Reset image size

Among the different strategies adopted to address the role of neutral impurities, the T-matrix approach to the Anderson problem is likely the best suited and most illuminating (Peres et al 2006, 2007a, Robinson et al 2008, Peres 2010, Wehling et al 2010).

The T-matrix has been introduced in section 2.5.2 and, for the present problem, has been given explicitly in terms of the unperturbed Green function G⁰ and the scattering potential $V$, irrespectively of its strength. According to scattering theory it determines the scattering amplitude and the scattering cross section (Taylor 1969). In particular, for Dirac fermions in 2D the differential cross-section can be given as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{{\rm d}\sigma ({\bf k,k'})}{{\rm d}\Omega} = \frac{A^2}{2\pi}\frac{k_{\rm F}}{\hbar^2 v_{\rm F}^2}|T_{{\bf k,k'}}|^2\equiv \frac{A^2}{2\pi}\frac{k_{\rm F}}{\hbar^2 v_{\rm F}^2}\frac{|t(\epsilon)|^2}{N_{\rm C}^2} \nonumber \end{align} \tag{ 7 }$$

where A is the area of the sample, with N_C carbon atoms, and T-matrix elements are taken between ${\bf k, k'}$ states of the unperturbed lattice, normalized on such area. This gives the total cross section—which is also the transport cross-section in this case—in the form

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \sigma = 2\frac{A_{\rm C}^2}{\hbar^2 v_{\rm F}^2}k_{\rm F}|t(\epsilon)|^2 = 2\pi \rho^0(\epsilon)\frac{A_{\rm C}}{\hbar v_{\rm F}} |t(\epsilon)|^2 \nonumber \end{align} \tag{ 8 }$$

where the factor of two comes from the final-state-sum over the two valleys, and $A_{\rm C}=A_c/2$ is the area per C atom. This scattering cross section can be used to estimate the scattering rate when n_i ($\nu_i$) impurities per unit area (per C atom) are present

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \tau^{-1} = v_{\rm F} n_i \sigma = \frac{2\pi}{\hbar} \nu_i |t(\epsilon)|^2 \rho^0(\epsilon) = -\frac{2}{\hbar}\nu_i \Im(t(\epsilon)) \nonumber \end{align} \tag{ 9 }$$

or, equivalently, the elastic mean free path $l_e=v_{\rm F} \tau$. The conductivity follows from the 2D Drude–Boltzmann result for isotropic dispersion¹⁴

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \sigma_{\rm DC}=\frac{2e^2}{h}k_{\rm F} l_e \nonumber \end{align} \tag{ 10 }$$

and reads as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \sigma_{\rm DC}=\frac{e^2}{2\pi^2}\frac{k_{\rm F} v_{\rm F}}{\nu_i |t(\epsilon)|^2 \rho^0(\epsilon)}. \nonumber \end{align} \tag{ 11 }$$

In the unitary limit ($V$ large, corresponding to the single vacancy case) the conductivity becomes

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \sigma_{\rm DC} = \frac{e^2}{\pi h}\frac{\nu_e}{\nu_i}\ln^2\left(\nu_e\right) \nonumber \end{align} \tag{ 12 }$$

(where $\newcommand{\e}{{\rm e}} \nu_e=\epsilon^2/\epsilon_c^2$ is the number of charge carriers per C atom) and displays a sublinear dependence on the carrier density. If the full structure of T is kept, on the other hand, the result is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \sigma_{\rm DC}=\frac{e^{2}}{\pi h}\frac{1}{\nu_{i}}\left[\left(\eta-{\rm sign}(\epsilon)\sqrt{\nu_{e}}\ln\left(\nu_{e}\right)\right)^{2}+\pi^2\nu_{e}\right] \nonumber \end{align} \tag{ 13 }$$

where $\newcommand{\e}{{\rm e}} \eta=\epsilon_c (\epsilon-\epsilon_{\rm ad})/W^2$ and the conductivity is no longer e–h symmetric, rather it presents smaller values for electrons (holes) when $\newcommand{\e}{{\rm e}} \epsilon_{\rm ad} > 0$ ($\newcommand{\e}{{\rm e}} \epsilon_{\rm ad} < 0$).

These expressions for the conductivity evidence the contribution of resonant states at or near the Dirac point. Their practical application only requires the tight-binding parameters of the Anderson problem and these are typically obtained by fitting DFT results, although some care is needed to obtain physically sound results (Robinson et al 2008, Wehling et al 2009). The ensuing resistivity $\rho_{\rm DC}=\sigma_{\rm DC}^{-1}$ provides the resonant scatterer contribution to the total resistivity, and should be considered along with other sources of disorder and electron–phonon scattering. The latter contribution has been found to be rather small by first principles calculations (Shao et al 2013), which estimated the room temperature intrinsic mobility to be as large as 3.3  ×  10⁵ cm² V⁻¹ s⁻¹ for both electrons and holes. Hence, with the exception of very cleaned samples such contribution can be safely neglected.

We focus here on an alternative approach designed to obtain the scattering cross section and the conductivity for neutral defects on graphene directly from first principle calculations (Achilli and Martinazzo 2018). This has the advantage of not requiring any parametrization, and to account self-consistently for any charge re-distribution that can occur around the defect. The method is similar in spirit to that introduced by Markussen et al to describe transport in silicon nanowires (Markussen et al 2007), and uses the non-equilibrium Green's function method (in conjunction with a Kohn–Sham Hamiltonian) to compute the Landauer's conductance of a scattering region containing one single defect. The reflection probability for traversing such region depends on both its width W and its length L, but while it saturates quickly when increasing L it is expected to decay as $\sim\sigma/W$ for increasing width. Hence, by performing a series of transport calculations for different dimensions of the scattering region—and eventually averaging over different configurations of the defect when necessary—one can obtain the desired scattering cross-section. The approach is 'semiclassical' because the very use of a cross-section implies that phase coherence is lost between a collision and the next.

To obtain the specific relation between the Landauer's transmission and the cross-section one may argue as follows. Let W_d and L_d be the transverse and longitudinal dimensions of the graphene device we are interesting in, and let n_w be the concentration of adatoms in such sample, $n_w=w^{-2}$, where w is average distance between them. Suppose that the defects are arranged on an array $W_d/w \times L_d/w$ of macroscopically small but microscopically large scattering regions, each containing one defect only. The sample conductance is that of W_d/w equal wires connected in parallel, each of which is a series of L_d/w equal resistors, namely

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\text G}_d = \frac{W_d}{w} \times \frac{w}{L_d} \times {\text R}_i^{-1} \nonumber \end{align} \tag{ 14 }$$

where ${\text R}_i$, the resistance of each circuit element, can be identified with the intrinsic resistance in the Landauer approach

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\text R}_i=\frac{h}{2e^2M}\left(\frac{1}{T}-1\right)=\frac{h}{2e^2M}\frac{R}{T} \nonumber \end{align} \tag{ 15 }$$

R and T being the average reflection and transmission coefficient per mode, and M the number of conduction channels available in the width w. It follows

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \sigma_{\rm DC}=\frac{2e^2}{h}M\frac{T}{R}\approx\frac{2e^2}{h}\frac{k_{\rm F} w}{\pi}\frac{T}{R} \nonumber \end{align} \tag{ 16 }$$

where k_F is the Fermi's momentum. Then, on comparing with the Drude–Boltzmann result $\sigma_{\rm DC}=\frac{2e^2}{h} k_{\rm F} l_e$, one obtains the elastic mean free path $l_e=\frac{w}{\pi} \frac{T}{R}$ and the transport cross-section

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \sigma = \frac{1}{n_w l_e} = \pi w \frac{R}{T}\label{eq: semicl conductivity}. \nonumber \end{align} \tag{ 17 }$$

In practice, equation (17) may be evaluated with first principles means for a number of (microscopic) graphene channels of increasing width $W \ll w$, and its limiting value inferred from the behavior of the sequence. The results of such calculations for the case of a H atom adsorbed on graphene are reported in figure 8, where the ab initio cross-sections are compared with those obtained with the T-matrix approach described above (Achilli and Martinazzo 2018). This figure shows that the ab initio cross section describes a much broader resonance than that found for the Anderson problem. The first principles cross section is larger than the model one in the entire energy range, with the exception of a narrow energy region centered at the position of the resonance in the model. The origin of such difference most likely lies in the fact that the first principles calculations account for the appropriate charge redistribution and, more importantly, for the structural changes that occur upon adsorption. In other words, the cross-section reflects both the changes in the electronics and the extended surface puckering (figure 1), a feature that the simple model of section 2.5 cannot describe.

Zoom In Zoom Out Reset image size

Once the cross section has been determined as described above the ab initio semiclassical conductivity for the given adatom concentration n_i can be calculated from the Drude–Boltzmann equation with the computed elastic mean free path

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \sigma_{\rm DC} = \frac{2e^2}{h}\frac{k_{\rm F}}{n_i\sigma}. \label{eq:conduct} \nonumber \end{align} \tag{ 18 }$$

In figure 9 the conductivity of two different graphene samples, one suspended (left panel, (Bolotin et al 2008)) and one supported on SiO₂ (right panel, (Novoselov et al 2005a)) have been reported along with the first principles semiclassical results for H. The experimental data in the left panel are well described by equation (18) using the computed first principles cross-section and assuming a concentration of 4.5  ×  10⁻⁶/C atom (1.7  ×  10¹⁰ cm⁻²) neutral scatterers. On the contrary, the data in the right panel required both neutral and charged scatterers, at concentrations of 3.0  ×  10⁻⁴/C atom (1.14  ×  10¹² cm⁻²) and 3.5  ×  10⁻⁴/C atom (1.33  ×  10¹² cm⁻²), respectively. Charged impurities are needed to improve the description for electrons ($V_g > 0$) while keeping that for holes at a reasonable level, and they were accounted for following equation (25) in Peres (2010). The presence of Coulomb scatterers is reasonable, since the sample of figure 9 (right) lies on a support that can trap charges, differently from the sample of the left panel, for which the presence of charges is hardly justifiable.

Zoom In Zoom Out Reset image size

The semiclassical approach for the calculation of conductivity is not suitable to describe the conductivity near the Dirac point where quantum interference effects become important (Adam et al 2009). In this energy region both experiments and theoretical predictions claim that graphene is characterized by a non-vanishing conductivity, despite the vanishing density of states (Peres 2010, Das Sarma et al 2011, Katsnelson 2012). While calculations for ideal graphene predict a value $\sigma_{\rm min}$ equal to $4e^2/\pi h$ the theoretical predictions for disordered graphene depends on the method adopted (Ziegler 1998, Katsnelson 2006, Tworzydło et al 2006, Ziegler 2006, Cserti 2007, Ziegler 2007) and the measured values are about two times larger (2e²/h) (Novoselov et al 2005a, Zhang et al 2005). This discrepancy stimulated an intense debate on the origin of the quasi-universal value of the conductivity minimum, ultimately leading to the conclusion that its value is related to the degree and nature of disorder in the carbon sheet (Suzuura and Ando 2002, Peres et al 2006, Schuessler et al 2009, Sajjad et al 2015). In particular, the minimal conductivity increases with disorder (Bardarson et al 2007, Titov 2007) and remains essentially constant when the temperature T_s is lowered by several orders of magnitude down to 30 mK (Tan et al 2007), where Anderson localization is expected (Das Sarma et al (2011) and reference therein). The absence of localization in graphene in the intermediate disorder regime indicates that the dominant disorder is either of long range character (Nomura and MacDonald 2007) and does not mix valley (Ostrovsky et al 2007) or preserves a chiral symmetry of the Hamiltonian (Ostrovsky et al 2010).

Hydrogen atoms, and resonant scatteres in general, have been envisioned as possible contributors to the minimum of conductivity (Stauber et al 2008) via the midgap states generated near the Dirac cones (Yuan et al 2010). These zero energy modes, that are proper of chiral symmetric graphene, are unaffected by Anderson localization and inelastic scattering effects, thereby giving rise to a robust metallic state characterized by a value of the conductivity independent of the impurity concentration (Ferreira and Mucciolo 2015), at least for ideal graphene with exact e–h symmetry.

The approach presented in the section 2.6.2, by addressing one (or a few) defects at a time, cannot describe localization effects, but intrinsically accounts for the quantum behavior of the electron dynamics in a complete ab initio way. Here, without attempting to justify the value of the minimum of conductivity that was the subject of several investigations (Peres et al 2006, Mucciolo and Lewenkopf 2010, Das Sarma et al 2011), we mention that our calculations do show a universal behavior at the Dirac point, i.e. one that depends on neither the dimensions of the model device used nor the nature of the impurity (being it a hydrogen or a fluorine atom, or a carbon atom vacancy). We shall address this issue in a future publication (Achilli and Martinazzo 2018).

The rationale for the preferential sticking mechanism is best given in terms of the RVB model introduced in section 2.5 (Casolo et al 2009a). Assuming that a first H atom adsorbs on an A-type site, the unpaired electron left on the B sublattice may easily (i.e. with a small or even vanishing activation barrier) couple to the electron of an incoming H to form a CH bond. Hence, formation of a 'AB dimer' is an easy process, and a singlet ground-state results in which the aromaticity of the substrate is partially restored. Conversely, if secondary adsorbtion occurred on the same sublattice to form an 'A₂ dimer', the incoming H atom would not take advantage of the available unpaired electron density (spin density), and adsorption would be as difficult as the first one. This is shown in figure 11 that reports the results of DFT calculations on a number of dimers as a function of the site-integrated magnetization, which is a measure of the average number of unpaired electrons available and roughly equals $|\psi({\bf r})|^2$ for the midgap state wavefunction at that site. Binding (barrier) energies are seen to increase (decrease) linearly by increasing the local magnetization, thereby showing how the spatial distribution of the midgap state determines the thermodynamic stability and the chemical reactivity. The result is a genuine electronic effect (in principle tunable by charge doping (Huang et al 2011)): substrate relaxation effects, though substantial, are site-independent for all but the so-called ortho dimer (with two Hs on neighboring sites), since surface puckering upon adsorption involves to a first approximation nearest-neighboring C atoms only. Recently, a detailed look at the relaxation effects has been provided by Zhang et al (2018) who used DFT calculations on large supercells to investigate the A–A and A–B interactions at long range. It has been found that the deformation potential makes the H–H interactions always attractive, though the A–B ones turn out to be much larger than the A–A one.

Zoom In Zoom Out Reset image size

The formation of dimers raises the issue of the competition between the ER reaction and secondary adsorption. Modeling the dynamics is quite challenging because it requires the development of a fully corrugated potential model that accounts for H₂ formation as well as for sticking at, at least, the ortho and para positions. The problem was solved by Casolo et al (2013) and (2016) in a classical though ab initio setting, i.e. with the help of AIMD. The advantage of the technique is that it avoids the need of computing and fitting a complicated potential, energies and forces being computed as required by the dynamics itself. The results of Casolo et al (2013) and (2016) showed the expected preferential formation of 'balanced' dimers, with the para dimers being the most abundant ones (up to one-half the total dimer fraction). Formation of these dimers features a dynamical threshold of  ∼20 meV (corresponding to a beam temperature $T_{\rm g}\sim 300$ K), below which only recombination is possible. This agrees well with the findings of Areou et al, who observed complete refreshment of a pre-covered surface upon exposing it to a low energy H atom beam (Aréou et al 2011).

Overall, ER abstraction dominates over dimer formation for collision energies smaller than  ∼0.2 eV. Nevertheless, accounting for the effects that formation of dimers may have in H₂ recombination was found to be crucial to improve the agreement between theory and experiments on the size of the reaction cross-section (Zecho et al 2002b). The cross-section obtained by Casolo et al of 3.1  ±  0.2 $\mathring{\rm A}^2$ at  ∼0.20 eV (corrected for the spin-statistical factor) agreed well with the value of  ∼3.5 $\mathring{\rm A}^2$ measured by Zecho et al at saturation for the same energy (that is, subtracting off the HA contributions, see section 2.4).

Low energy orbitals show a marked tendency to localize at the edge sites, irrespective of whether or not they are regularly shaped. This is a consequence of e–h symmetry, i.e. of the symmetry of the tight-binding Hamiltonian with nearest-neighbors interactions only, $H=H_{AB}+H_{BA}$ (see section 3.4). To show this, we perform a lattice 'renormalization' (Naumis 2007, Martinazzo et al 2010) and focus on one sublattice only (say $\mathcal{A}$) and on the 'renormalized' Hamiltonian $\tilde{H}_{AA}=H_{AB}H_{BA}$. Indeed, since H only allows transitions from the A to the B subspaces (H_BA) and vice versa (H_AB) it is sufficient to consider the problem in the A space only with the Hamiltonian¹⁷ $\tilde{H}_{AA}$.

The spectral properties of $\tilde{H}$ are closely related to those of the original Hamiltonian H. For any non-zero eigenvalue $\newcommand{\e}{{\rm e}} \tilde{\epsilon}_{i}$ and eigenvector $\renewcommand{\ket}[1]{|#1\rangle} \ket{\psi_{A, i}}$ of this Hamiltonian there exist two solutions of the original problem with eigenvalues $\newcommand{\e}{{\rm e}} \epsilon_{i}=\pm\sqrt{\tilde{\epsilon}_{i}}$ and eigenvectors $\renewcommand{\ket}[1]{|#1\rangle} \ket{\psi_{A, i}}\pm\ket{\psi_{B, i}}$, where $\renewcommand{\ket}[1]{|#1\rangle} \ket{\psi_{B, i}}$ is defined to be $\renewcommand{\ket}[1]{|#1\rangle} \newcommand{\e}{{\rm e}} \ket{\psi_{B, i}}=\tilde{\epsilon}_{i}^{-1/2}H_{BA}\ket{\psi_{A, i}}$ (if $\newcommand{\e}{{\rm e}} \tilde{\epsilon}_{i}=0$, $\renewcommand{\ket}[1]{|#1\rangle} \ket{\psi_{A, i}}$ is already a H eigenvector). The converse is also true, namely from any eigenvector $\renewcommand{\ket}[1]{|#1\rangle} \ket{\psi_{i}}$ the two projections $\renewcommand{\ket}[1]{|#1\rangle} \ket{\psi_{A, i}}$ and $\renewcommand{\ket}[1]{|#1\rangle} \ket{\psi_{B, i}}$ onto the A and B subspaces satisfy $\renewcommand{\ket}[1]{|#1\rangle} \newcommand{\e}{{\rm e}} H_{BA}\ket{\psi_{A, i}}=\epsilon_{i}\ket{\psi_{B, i}}$ and $\renewcommand{\ket}[1]{|#1\rangle} \newcommand{\e}{{\rm e}} \tilde{H}_{AA}\ket{\psi_{A, i}}=\epsilon_{i}^{2}\ket{\psi_{A, i}}$; that is, in studying $\tilde{H}_{AA}$ one only misses possible zero eigenstates in the B subspace, that can be easily detected by defining A to be the majority species and comparing the number of zero A eigenstates with the sublattices imbalance. Importantly, the above lattice renormalization maps the low-energy side of the excitation spectrum of the original problem into the low-energy sector of $\tilde{H}_{AA}$: the ground-state of the renormalized system corresponds to the highest occupied/lowest unoccupied molecular orbital (HOMO/LUMO) pair of the original system—the so called frontier orbitals. Frontier orbitals are key quantities because, from a chemical perspective, they govern the system chemical reactivity to a large extent (Fukui et al 1954), similarly to what midgap states do when unpaired electrons are present. Hence, the ground-state of the renormalized problem is a powerful tool to investigate the reactivity, a quite unique situation.

The renormalized lattice of a honeycomb lattice (see section 2.5.2, with NN hoppings only) is a hexagonal (commonly referred to as 'triangular') lattice, that is of course the sublattice $\mathcal{A}$ of the original system, with different parameters. The renormalized hopping is just t² (assuming t_ij  =  t for simplicity), while the on-site energies are $t^{2}Z_{i}$, where Z_i is the coordination number of the i^th A site in the original lattice, a sort of 'π-coordination number'. In the bulk of graphene Z  =  3 but at an edge such number is smaller, typically Z  =  2, though Z  =  1 is also possible¹⁸. Of course, the ground-state of the renormalized lattice tends to localize on the sites with the lowest on-site energy; hence, frontier orbitals of the original π system must show a certain degree of edge localization. Furthermore, among the edge sites, the ground-state of the fictitious problem localizes mainly on the sites with the largest number of similar neighbors in the renormalized lattice—next-to-nearest neighbors in the original lattice—since this situation maximizes the degree of hybridization. Hence the number (ξ) of next-nearest-neighbors provides further insights into the nature of the edge state (Bonfanti et al 2011b). We named it the hypercoodination number since, unfortunately, the more appropriate term hyperconjugation was already in use in the chemical literature, with a rather different meaning. Its usefulness becomes soon evident when considering the textbook cases of armchair and zig-zag edges. The 'true' edge sites are two-fold π coordinated in both cases (Z  =  2) but the hypercoordination is larger for the zig-zag ($\xi=2$) than for the armchair ($\xi=1$) one, and this suggests a more marked edge localization (and reactivity) for the first.

The combination of edge localization and sublattice imbalance (with its own sublattice localization) provides a set of useful rules that helps predicting graphene reactivity. Two different situations can be envisaged that alternate to each other in a sequential adsorption process at the edges (Bonfanti et al 2011b, Jensen et al 2018).

In the absence of unpaired electrons reactivity is largest at the edge sites with lowest Z and largest ξ. The prototypical cases is H atom adsorption at a perfect (H-terminated) zig-zag edge, where the binding energies is $E_b \sim 2.8$ eV, i.e. much larger than for a 'bulk' site ($E_b\sim 0.9$ eV, see section 2.1). DFT (and higher theory level) calculations on a number of graphene clusters of different size and shape showed that the H binding energy at the Z  =  2 sites of irregularly shaped edges fall all between the above mentioned limits, with a clear preference for sites with largest ξ (Bonfanti et al 2011b). Because of its nature the ξ number may fail when comparing sites in different systems, but it is rather robust when used within the same structure. In addition, reliable upper bounds for the binding energy exist: the one given above for perfect zig-zag edges is the largest (useful for $\xi=2$), smaller upper bounds are  ∼1.7 eV for $\xi=1$ (the binding energy for the perfect armchair edge) and  ∼1.3 eV for $\xi=0$. For 'bulk'-like graphenic sites (Z  =  3), on the other hand, ∼1.0 eV can be considered a reliable bound.

When the adsorbed hydrogen atoms leave a sublattice imbalance on the surface, one or more unpaired electrons (midgap states) appear. Similarly to those discussed in section 5, such states localize on the majority sublattice but they are now shaped by the presence of the edge, showing a preference for small Z and large ξ sites, for the same reasons given above. Binding hydrogen atoms on such sites is rather easy and binding energies can be rather large, up to  ∼3.0 eV for Z  =  2, ∼3.8 eV for Z  =  1 and  ∼4.0 eV for Z  =  0 (the binding energy for forming the para dimer provides the right figure for the case Z  =  3, $E_b\sim2.0$ eV).

As a consequence, hydrogen atoms—when employed under mild reaction conditions that prevent bulk adsorption—should start binding at the edges and 'corrupt' the graphene sheet by reducing the spatial extension of its π cloud. In this way, as hydrogenation proceeds carbon atoms that were initially in the bulk (hence relatively inert) can find themselves at the 'moving' edges of the π cloud (Z  =  2) or be even isolated from it ($Z=0, 1$), thereby becoming rather prone to hydrogenation. That reaction can indeed proceed starting from the edges has been shown with Birch-type 'wet' hydrogenation (Zhang et al 2016) on mechanically exfoliated graphene (figure 12). Birch reduction is a hydrogenation reaction that has been long used to partially hydrogenate aromatic rings, e.g to convert aromatic rings into 1,4-cyclohexadienes. It differs from atomic hydrogen exposure¹⁹ but it involves yet the spatial properties of the frontier orbitals discussed above, since they are required to accommodate an electron from the reaction medium. Interestingly, only bilayer and few-layer graphene have shown edge-hydrogenation, the more reactive single layer graphene was seen to undergo uniform reaction.

Zoom In Zoom Out Reset image size

Compelling evidence for silicene, the two-dimensional allotrope of Si, was first reported in 2012 (Vogt et al 2012). This material has attracted much attention from the community since the beginning, mainly because of the expectations of being easily integrated into the existing Si-based technology, partly met with the realization of the first silicene field effect transistors (Tao et al 2015). Silicene has a band structure very similar to graphene, with two Dirac cones at the corners of its hexagonal Brilluoin zone, and a very similar Fermi velocity. At a closer look, however, there exist marked differences. Silicene is not flat, rather presents a buckled structure that makes the sublattices distinguishable by e.g. application of an electric field. Hence, the Si atoms have some sp³ character on their own. It has further (much) weaker π (and σ) bonds, that implies that its lattice is softer than graphene. Thus, it is instructive to consider the hydrogenation process of silicene, and compare it with the situation in graphene (Pizzochero et al 2016).

Hydrogen adsorbs on silicene much like it does on graphene: the hydrogen atom binds on top a Si atom and triggers a (further) outward motion of the lattice atom. The energetics, though, are completely different: the binding energy is $E_b\sim2.2$ eV, more than 1 eV larger than in graphene. This is the combined effect of structural and electronic factors. Since the lattice is softer, the reconstruction energy (the energy going into the substrate to create a binding Si with a full sp³ character) is only  ∼0.1 eV, one order of magnitude less than in graphene, and this energy gained is made available as binding energy. Secondly, weaker π bonds mean that less energy is needed to break them in order to form a bond with the adatom. Noteworthy, hydrogen adsorption on silicene turns out to be barrierless, hence much easier than in graphene. This has important practical consequences: silicene has similar transport properties of graphene, on account of their similar band-structure, but a reduced mobility is expected for it on the basis of its easier chemistry, i.e. the larger amount of adatoms that can be gathered in the fabrication process. In addition, dimer and clusters are thermodynamically favored as in graphene but the absence of a sticking barrier makes their presence relevant only under equilibrium conditions, when the Si sheet is left to equilibrate in a hydrogen atmosphere. The adsorption rate has the same high value for any site—the one expected for barrierless sticking—hence, under typical non-equilibrium conditions, hydrogenation of silicene is a rather random process, without any clustering tendency of the adatoms, in sharp contrast to graphene.

Single–layer hexagonal boron nitride (h-BN) is an another 2D material that has attracted much attention after the discovery of graphene. h-BN (aka graphitic BN) is the most stable polymorph of boron nitride, and is made of weakly (van der Waals) bonded planes, each composed by boron and nitrogen atoms arranged in a honeycomb lattice. Single–layer h-BN is isoelectronic to graphene and, similarly to it, has been obtained by mechanical exfoliation (Novoselov et al 2005b, Pacilé et al 2008, Lee et al 2010), epitaxial growth (Oshima and Nagashima 1997), ultrasonication (Zhi et al 2009) and high-energy electron beam irradiation of BN particles (Jin et al 2009). Chemical synthetic methods have also been successfully applied to obtain samples of larger and larger sizes, and with improved quality (Shi et al 2010, Lee et al 2012, Kim et al 2013, Lu et al 2015). The striking difference of single-layer h-BN with respect to graphene is its insulating behavior (Britnell et al 2012), which is due to a large direct band-gap in its electronic structure (Watanabe et al 2004, Cassabois et al 2016). From the perspective of the low-energy physics, this is a consequence of the 'mass' term in the Dirac-like equation that arise when the two sublattice sites of the honeycomb network are occupied by inequivalent atoms.

Hydrogenation of h-BN has been studied long before the emergence of the interest in 2D materials. One of the most basic problems considered is whether the adsorption of H happens preferentially at one of the two atomic species composing the surface. In light of the contradictory theoretical results present in the literature (Koswattage et al (2011) and references therein), two recent experimental works addressed this issue. Koswattage et al (2011) examined a deuterated two-layer sample of h-BN with NEXAFS and XPS, and proved that binding occurs mostly at the B sites. This conclusion was further confirmed by Ohtomo et al, who carried out a more comprehensive experimental and theoretical study, and considered the hydrogenation of a h-BN layer grown epitaxially on Ni(1 1 1) (Ohtomo et al 2017). Apart from the selectivity towards the B sites, however, binding of H occurs with an electronic and structural reorganization similar to that found in graphene. Upon hydrogenation the B atom undergoes a similar re-hybridization, leading to a partial tetrahedralization of the lattice.

One further aspect that has been highlighted in epitaxial h-BN is the possibility of intercalating atomic hydrogen rather than adsorbing it on its surface. Brugger et al, for instance, monitored the Moiré pattern that h-BN forms on Rh(1 1 1) during exposure to hydrogen atoms and found that H intercalation causes an effective decoupling between the 2D layer and the substrate (Brugger et al 2010). Similar results have been obtained for h-BN on Ni(1 1 1), although in this case intercalation takes place at a higher level of H exposure and is less stable than hydrogenation (Späth et al 2017).

It is worth noticing that, similarly to graphene, hydrogen adsorption has been investigated in light of the band-gap modulations that would be needed for applications in electronic devices. Unlike graphene–related substrates, however, the goal in this case is the reduction of the band–gap, and the conversion of h-BN to a semi-conductor. Along this line, some interesting results have been obtained by Zhang and Feng, who observed a significant increase of conductivity after exposing rippled membranes of few–layer h-BN to a hydrogen plasma (Zhang and Feng 2012).

Among the emergent 2D materials, we also mention phosphorene, an allotrope of phosphorus that can be viewed as a single layer of P atoms, although arranged effectively in three–dimensions. Contrary to the other single–layer materials just examined, hydrogenation of phosphorene has been rarely considered. However, there exist a few theoretical investigations that examined adsorption of H as one of the possibilities for doping the substrate with heteroatoms (Boukhvalov et al 2015, Kulish et al 2015).

Single layer- (SL), bilayer- (BL) and few layer- (FL) graphene are known to have different electronic structures, that differ also from that of graphite (Castro Neto et al 2009). Nevertheless, it is often assumed that they are all very similar to each other from a chemical point of view, since the interaction between layer is so a small fraction of the H–graphene interaction that it hardly modifies the physics of the adsorption process. As mentioned in the introduction we too made such an assumption in this manuscript and often used graphite as a graphene knockoff, when data on graphene were unavailable. It is thus important to address this issue in detail, and highlight the differences between free-standing SL graphene and graphene interacting with one or more graphene layers.

Apart from the role of curvature discussed above, one key issue in this context is the strength of the interaction between the layers, i.e. the binding energy in BL and/or one of the several (slightly different) interlayer-interaction energies that can be defined in multilayer graphene/graphite. The binding energy in BLG is not experimentally known, but its value is bound from above by the interlayer cohesive energy of graphite—the energy needed to separate the layers, over the total number of C atoms—that was determined to be 52  ±  5 meV/atom (Zacharia et al 2004). Plane-wave vdW-inclusive DFT calculations indicate that the graphite cohesive energy is ca. 10% larger than the binding energy in BLG (see below), so additional layers have a minor importance. In the situation we are interested in, such interlayer interaction represents a restoring force opposing to the surface puckering that accompanies H adsorption, hence it is expected to reduce the hydrogen binding energy when going from SLG to BLG, FLG and graphite. Furthermore, depending on the layer stacking it may also break the sublattice symmetry. In the typical situation of an AB stacking, it determines a 'chemisorption sublattice imbalance' already in BLG, since it acts differently depending on whether the H atom adsorbs on the 'softer' β site (the sites without a C underneath) or on the 'stiffer' α site (the sites with a C underneath). This has led to the fascinating hypothesis that a transient sublattice imbalance (hence a transient magnetization) can be realized when thermally desorbing hydrogen atoms from BL- and ML- graphene (Moaied et al 2015).

DFT calculations with vdW-inclusive functionals, properly corrected for the BSSE if necessary, reproduce the available experimental data on the interlayer interactions rather well and can thus reliably predict the energetics for hydrogen adsorption (Bonfanti and Martinazzo 2018). Table 1 gives some significant energies (in meV per atom) in BLG and graphite, as obtained with two different vdW-inclusive density functionals, vdW-DF (Dion et al 2004, Kong et al 2009) and vdW-VV (Vydrov and Van Voorhis 2010). As can be seen from the table, the two functionals perform similarly well for the interlayer interaction, but the second turns out to be superior when investigating H adsorption (Bonfanti and Martinazzo 2018). The resulting chemisorption energy on SLG is larger than the values it takes on BLG, in agreement with expectations and experimental observations. The reduction of binding energy when passing from SLG to BLG compares favorably with the interlayer interaction energy (per atom), as it should since this is the energy lost by the binding carbon atom when it is pulled out of the surface in the adsorption process. The site-dependency of the chemisorption energy turns out to be smaller, of the order of some tens of meV.

It is worth stressing at this point that addressing vdW interactions is yet challenging nowadays from a theoretical point of view. It is true that recent, vdW-inclusive density functionals can achieve impressive performances, but their reliability is not yet general, and needs to be checked on a case-by-case basis. A common problem that is becoming more and more frequent, though, is the use of the vdW-inclusive functionals in conjunction with an atomic-orbital-based DFT approach, a rather dangerous mix that is largely overbinding because of the basis set superposition error. In the case of H on graphene, for instance, this overestimation of the vdW interactions may completely wash the sticking barrier (Moaied et al 2015, Brihuega and Yndurain 2018), clearly at odds with the many evidences brought above on its existence and its consequences.

On the experimental side, higher reactivity of SL towards hydrogenation (w.r.t. to ML graphene) has been reported by several authors, using different hydrogen sources and graphene supports (e.g. Ryu et al (2008), Elias et al (2009) and Zhang et al (2016)) though opposite findings have also reported (Luo et al 2009). This suggests that, in practice, the preference for SL graphene is probably subtle and can be reverted by changing the operation conditions (e.g. the H atom source).

Graphene has been grown on a huge number of metal substrates (Wintterlin and Bocquet 2009, Tetlow et al 2014, Dedkov and Voloshina 2015), either by surface segregation or deposition of hydrocarbons. Perfectly ordered (although often misaligned) epitaxial overlayers can be obtained on a number of hexagonally close-packed surfaces, where Moiré structures with large periodicities appear, as a consequence of a small lattice mismatch between graphene and the underlying metal. The metal–graphene spacings can vary in a wide range between 2.1 and 3.8 Å depending on the strength of the interaction between the carbon sheet and the substrate, that also reflects on the graphene electronic properties. In some systems the $\pi/\pi^*$ bands are significantly hybridized with the metal band structure. In other systems graphene is quasi free-standing and displays an almost intact electronic structure. Hydrogenation of the carbon sheet can then be performed in different chemical environments, depending on the specific metal susbtrate chosen. Reviewing the literature on such a topic is well beyond the aims of the present work. Rather, in the following, we use a few examples that illustrate the situations typically found in practice.

Depending on the extent of the graphene–metal interactions, we may distinguish four different 'classes'. In a first, graphene is and remains weakly bound to the substrate, even upon extensive hydrogenation. This is rather similar to hydrogenating free-standing graphene, with little (if any) effect of the underlying metal substrate. In a second class fall those combinations where the graphene overlayer is initially quasi free-standing, but shows an increased interaction with the metal substrate upon hydrogen adsorption. A third class is for those combinations where graphene is and remains strongly bound to the surface, while a fourth (apparently yet empty) is for graphene layers that are initially strongly bound to the metal substrate but progressively freed from it as hydrogenation proceeds. Clearly, the progress of the hydrogenation process, as well as the ordered structures that hydrogen can form, may vary largely depending on the graphene-substrate interactions, both 'early' and 'late'. The simplest thing to consider, for instance, is to what extent the induced magnetic moment of a H adatom gets quenched by the substrate, since this determines whether orienting effects innate in the carbon sheet are to be expected or it is the substrate that plays a primary role in the arrangement of the adatoms.

'Graphene on gold' belongs to the first class. Graphene binds strongly on Ni(1 1 1) but intercalation of Au progressively reduces its interaction. 1 ML of Au entirely decouples the carbon sheet from the substrate and makes it quasi free standing and undoped (Varykhalov et al 2008). Angle-resolved photoemission spectroscopy revealed indeed a gapless linear π-band dispersion near K and a Dirac point within 25 meV to the Fermi level. In addition, the Au layer underneath the carbon sheet provides a chemical inert buffer layer that prevents 'late' interactions with the susbtrate upon functionalization. This makes Gr/Au/Ni(1 1 1) an ideal system where investigating hydrogen reactivity (Haberer et al 2010, 2011a, 2011b). Several controlled hydrogenation experiments have been conducted on such system, and saturation values of 25% observed. Furthermore, kinetic experiments have been also performed with both hydrogen and deuterium, and deuteration has been found to proceed faster than hydrogenation, and to lead to substantially higher saturation coverages of 35% (Paris et al 2013).

Ir(1 1 1) is a substrate which couples weakly to the graphene sheet, as evidenced by the presence of characteristic Dirac cones in the quasi particle spectra (Pletikosić et al 2009). Graphene lays at $3.38\pm0.04$ Å above the metal (Busse et al 2011), and it is only weakly corrugated (Hämäläinen et al 2013), with a variation of its height of $0.47\pm0.05$ Å over the  ∼25 Å periodicity of the Moiré pattern that originates from the 10% lattice mismatch with Ir(1 1 1) (see figures 14(a) and (b)). Graphene reactivity though varies over the surface. In the so-called hcp and fcc areas of the Moiré, where the position of every second carbon atom coincides with the position of an Ir atom below (see figure 14(d)), the graphene lattice can distort upon hydrogen functionalization, and bind stronger to the metal. This is already clear with a single hydrogen atom (see figure 14(c)) and becomes more plausible when multiple hydrogenation is considered, since a configuration is possible in which every second C atom binds to the underlying Ir, while neighboring C atoms bind to H atoms on top. A selective functionalization of hcp areas by hot H atoms has been achieved and seen to give rise to highly ordered hydrogenated structures, including the opening of a gap in the electronic band structure (Balog et al 2010, 2013, Jørgensen et al 2016). Graphene/Ir(1 1 1) is interesting in many respects, as it has been recently shown that the carbon sheet can 'mediate' the catalytic activity of the substrate, and make dissociative adsorption of H₂ feasible when employing vibrationally excited molecules (Kyhl et al 2018). In this case, functionalization of the graphene surface occurs in a highly ordered manner and exhibits an avalanche effect where the first dissociative adsorption event decreases the barriers for subsequent dissociative adsorption.

Zoom In Zoom Out Reset image size

Finally, graphene/Ni(1 1 1) represents a case of strong interaction between the carbon sheet and the substrate. This and similar other systems have long been investigated in surface science, where they were found to exhibit a high crystalline quality and inertness in air (Gamo et al 1997). On Ni(1 1 1) a downshift of the Dirac point occurs that causes the hybridization of the graphene π bands with the Ni 3d states. As a consequence, the original electronic structure of graphene at low energy is fully destroyed. Despite this, graphene on Ni(1 1 1) is essentially flat, and lies at about 2.1 Å above the surface (Gamo et al 1997). It forms a 1  ×  1 overlayer with two inequivalent carbon atoms, one on top of a Ni atom and one on a hollow site²⁰. Similarly to the case discussed above, the structural changes upon functionalization favor C–H formation on the carbon atoms located on the hollow sites, since in this way the neighboring carbon atoms at top sites buckle downward and increase the strength of the carbon-substrate bonds (Andersen et al 2012). This structural change also affects the adsorption of the subsequent hydrogen atoms, and meta dimers become favored over 'traditional' AB structures. Increasing the coverage, the half-hydrogenation of graphene with full occupation of one sublattice only becomes possible. This is the so-called 'graphone' that has been recently synthetized on Ni(1 1 1) (Zhao et al 2015).

The presence of a metallic susbtrate beneath the carbon sheet allows modulation of the chemical properties, and selective functionalization (hydrogenation) of specific area of the graphene sheet becomes possible. The realization of such patterned hydrogenation is of great technological interest, especially in relation to the long-standing issue of a band gap-opening. That such modulation is indeed operative was shown by Balog et al (2010) on the graphene/Ir(1 1 1) system discussed in the previous section. The authors of Balog et al (2010) were able to hydrogenate the hcp and fcc regions of the Moiré and to open a sizable bandgap in graphene (⩾0.5 eV) by effectively confining the charge carriers in a 'antidot' structure. The approach was later refined (Balog et al 2013), up to the point where selective hydrogenation of the hcp regions only was achieved by controlling the surface temperature during the hydrogenation (Jørgensen et al 2016).

Graphene/Au/Ni too has been shown to give rise to well ordered structures and a gap opening  ⩾0.5 eV in graphene (Haberer et al 2010, 2011b). Furthermore, a crystalline compound with 4:1 stoichiometry ('C₄H') was apparently realized (Haberer et al 2011a) that is the hydrogen analogue of the single-sided fluorinated product obtained from graphene on Cu(1 1 1) (Robinson et al 2010). These structures are the simplest of a general class of graphene superlattices where a band gap opens because of the preservation of symmetry (Martinazzo et al 2010). Such superlattices have sizable gaps that scale optimally with the superlattice constant, and, in addition, preserve Dirac carriers at the gap edges.

We further notice that patterned adsorption can also be obtained via appropriate treatment of the graphene surface. For example ordered lines of hydrogen dimers have been obtained on graphite upon pre-covering the surface with cyanuric acid (Nilsson et al 2012). The highly stable ordered hydrogenated domains can be exploited to confine electrons similarly to nanoribbons and to realize chemical electron waveguides, as demonstrated by transport calculations (Achilli et al 2014).