## Electronic Transport in Graphene with Aggregated Hydrogen Adatoms

Fernando Gargiulo, 1 Gabriel Autès, 1 Naunidh Virk, 1 Stefan Barthel, 2, 3 Malte Rösner, 2, 3 Lisa R. M. Toller, 1 Tim O. Wehling, 2, 3 and Oleg V. Yazyev 1

1 Institute of Theoretical Physics, Ecole Polytechnique Fédérale de Lausanne (EPFL), CH-1015 Lausanne, Switzerland 2 Institut für Theoretische Physik, Universität Bremen, Otto-Hahn-Allee 1, D-28359 Bremen, Germany 3 Bremen Center for Computational Materials Science, Am Fallturm 1a, D-28359 Bremen, Germany (Dated: August 4, 2014)

Hydrogen adatoms and other species covalently bound to graphene act as resonant scattering centers affecting the electronic transport properties and inducing Anderson localization. We show that attractive interactions between adatoms on graphene and their diffusion mobility strongly modify the spatial distribution, thus fully eliminating isolated adatoms and increasing the population of larger size adatom aggregates. Our scaling analysis shows that such aggregation of adatoms increases conductance by up to several orders of magnitude and results in significant extension of the Anderson localization length in the strong localization regime. We introduce a simple definition of the effective adatom concentration x glyph[star] , which describes the transport properties of both random and correlated distributions of hydrogen adatoms on graphene across a broad range of concentrations.

PACS numbers: 72.80.Vp, 73.20.Hb, 71.23.An

Graphene has unveiled a plethora of unconventional transport phenomena [1-3], such as the universal minimal conductivity [4], Klein tunneling [5] and the anomalous quantum Hall effect [6, 7]. On the applied side, graphene is interesting because of its exceptionally high charge-carrier mobility, which is typically limited by the presence of various types of disorder. Resonant scattering impurities, such as chemical functionalization defects [8] and dislocations [9], show the most pronounced effects on charge-carrier transport in graphene. Hydrogen adatoms represent a prototypical resonant scattering impurity, which can be experimentally introduced in a controlled fashion [10] and allows for a simple theoretical description [11]. A hydrogen adatom covalently binds to a single carbon atom of graphene resulting in rehybridization into the sp 3 state, thus effectively removing that site from the honeycomb network of p z orbitals. This gives rise to a zero-energy state localized around the defect, and results in the resonant scattering of charge carriers.

At a fundamental level, the classical scaling theory of Anderson transition predicts complete localization of the electronic spectrum in two dimensions (2D), regardless of the amount of disorder [12]. For hydrogenated graphene, a model based on massless Dirac fermions with δ -function-like potentials confirms this prediction of the unitary class, though in 2D systems localization lengths can be strongly energy-dependent and, eventually, very large [13]. However, no unanimous consensus has been reached since experiments on hydrogenated graphene point towards metal-insulator transition, theoretically justified by the presence of electron-hole puddles (2D percolation class) [14-17].

Early works treating finite concentrations of resonant impurities in graphene assumed that the total scattering cross-section deviates little from the incoherent addition of the individual cross-sections, for example in the Boltzmann equation framework [18]. This picture is valid for low defect concentrations, low charge-carrier densities and random adatom distributions. A better description requires including the effect of coherent superposition of wavefunctions scattered by distinct adatoms [8, 19-21]. This is particularly important when impurities are in proximity to each other, with the limiting case being the formation of compact clusters in which hydrogen adatoms populate neighboring carbon atoms [22-24]. Indeed, the overall short-range attractive interaction between individual hydrogen adatoms on graphene [25, 26], combined with their relatively high diffusion mobility at room temperature [27-29], suggests a high degree of spatial correlation between adatoms.

In this Letter, we address the effects of spatial correlation of resonant impurities on electronic transport in graphene. The equilibrium configurations of hydrogen adatoms on graphene, obtained by means of MonteCarlo simulations, show a strong tendency towards aggregation into small clusters essentially eliminating isolated adatoms. Electronic transport properties investigated using the Landauer-Büttiker approach complemented with kernel polynomial method calculations show that aggregation dramatically increases both the conductivity and the localization length. We propose a unified framework to account for the effects of spatial correlation of resonant scattering centers on electronic transport in graphene.

Upon adsorption, a hydrogen adatom covalently binds to a single carbon atom of graphene changing its hybridization state to sp 3 and its coordination sphere to tetrahedral as shown in Fig. 1(a) [30]. The covalent binding of a second adatom to the nearest neighbor carbon atom partially releases the elastic energy due to the change of coordination sphere, thus resulting in effective attractive interaction between adatoms [25, 26]. This suggests that the interaction between adatoms can be

Figure 1: (a) Atomic structure of a hydrogen adatom covalently bound to graphene. (b) Structures of small clusters of hydrogen adatoms used for fitting the pair potential of Eq. (1) (c) Predicted energy ˜ E of aggregation of hydrogen adatoms as a function of aggregation energy E DFT calculated from first principles for the set of adatom clusters shown in panel (b).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0245v1-figure-d5517b1f-1eb4-4574-989f-335c631d13c1.jpg)

accurately described using a short-ranged pair potential. In our study, we expand an Ising-like interaction energy E up to the second nearest-neighbor term

$$E = \gamma _ { 1 } \sum _ { \langle i, j \rangle } s _ { i } s _ { j } + \gamma _ { 2 } \sum _ { \langle \langle i, j \rangle \rangle } s _ { i } s _ { j }, \quad \ \ ( 1 ) \quad \ s u b _ { \text{tain} }$$

where γ 1 and γ 2 are the corresponding first and second nearest-neighbor parameters. Here, s i = 1 if a carbon atom i is populated by an adatom, otherwise s i = 0 . Parameters γ 1 and γ 2 are obtained by fitting the interaction energies of adatoms calculated from first principles [31] for a set of small adatom aggregates shown in Fig. 1(b). Under the assumption of single-side functionalized graphene, the obtained parameters γ 1 = -1 182 . eV and γ 2 = 0 484 . eV [32] signify a considerable firstnearest-neighbor attraction alongside a weaker secondnearest-neighbor repulsion. The excellent agreement between the interaction energy ˜ E , estimated using the fitted γ 1 and γ 2 , and the first-principles values E DFT , confirms the applicability of the short-range pair potential form (1) for describing small clusters [Fig. 1(c)].

In order to assess the effect of interaction between adatoms on their spatial distribution and transport properties we perform Monte-Carlo simulations using the introduced pair potential (1). The simulations employ a Monte-Carlo move based on the displacement of a randomly selected adatom to a random unoccupied carbon atom in combination with the Metropolis acceptance criterion [33]. We considered models containing up to N C = 10 6 carbon atoms ( 165 × 165 nm ) and adatom con2 centrations x = N /N H C ranging from 0.1% to 10%. All simulations have been performed at T = 300 K . A representative configuration of randomly distributed adatoms [Fig. 2(a)] is compared with a configuration obtained from a Monte-Carlo simulation [Fig. 2(b)] at x = 5% .

Further details are revealed by comparing the cluster size distributions P n ( ) [Fig. 2(c)], with adatoms populating neighboring carbon atoms being assigned to the same cluster. In the case of a random distribution most adatoms are isolated ( n = 1 ), while the occurrence of clusters ( n &gt; 1 ) is merely a probabilistic effect. In contrast, no isolated adatoms are found in the presence of interactions, with the most abundant species being adatom dimers ( n = 2 ). The size distribution for the correlated case shows a longer tail with a significant probability of observing up to n = 6 clusters. The dependence of P n ( ) on adatom concentration x is relatively weak [Fig. 2(d)].

We now focus on electronic and transport properties calculated using the nearest-neighbor tight-binding Hamiltonian for p z orbitals

$$H = - t \sum _ { \langle i, j \rangle } [ c _ { i } ^ { \dagger } c _ { j } + \text{h.c.} ]$$

with the hopping integral constant t = 2 7 . eV [2]. An adsorbed hydrogen atom is modeled by excluding the p z orbital of the carbon atom to which it is bound as a consequence of sp 3 hybridization, making it similar to a vacancy defect [27, 30, 34]. We stress that adatoms do not introduce coupling between the sites owing to the same sublattice of the bipartite lattice of graphene, thus maintaining electron-hole symmetry of the electronic spectrum. Figures 3(a,b) compare the density of states (DOS) of graphene with random and correlated distributions of adatoms at different concentrations. In the case of a ran-

Figure 2: (Color online) Representative configurations of (a) randomly distributed and (b) correlated hydrogen adatoms on graphene at x = 5% concentration. (c) Comparison of random and correlated adatoms cluster size distributions P x ( ) at x = 5% concentration. (d) Cluster size distributions P x ( ) of correlated hydrogen adatoms at different concentrations. All correlated configurations are obtained by means of MonteCarlo simulations carried out at T = 300 K.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0245v1-figure-47e07739-4c72-4bac-b310-1de6a5d4477a.jpg)

Figure 3: (Color online) Density of states of graphene in the presence of (a) randomly distributed and (b) correlated hydrogen adatoms at different concentrations. Local density of states (LDOS) on carbon atoms in the vicinity of (c) an isolated hydrogen adatom and (d) a dimer of hydrogen adatoms. In panels (c,d) LDOS referred to as 'analyt.' have been obtained using the analytical Green's function calculations [36].

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0245v1-figure-4bfa5eb9-f166-4c46-8e54-59898d50ee7d.jpg)

dom distribution one observes a strong peak at E = 0 due to the resonant modes originating from isolated adatoms [19, 34]. The corresponding wave function is localized on the sublattice opposite to that of the carbon atom binding the adatom and decays from the defect position [35]. At high adatom concentrations, x &gt; 1% , the E = 0 peak is accompanied by flat density regions at higher energies with a noticeable overall renormalization of the DOS, in agreement with previous calculations [34]. In comparison, the DOS calculated for the correlated impurity configurations shows a less intense peak at E = 0 and an increased weight for -0 9 . t &lt; E &lt; 0 9 . t that is more evident at higher concentrations. This change is a direct consequence of different cluster size distributions. The dominant cluster type in the case of the correlated adatom distribution is the dimer ( n = 2 ), which is known to be non-resonant [11, 24, 36], meaning that no localized states emerge at any energy. The local density of states (LDOS) calculated on neighboring atoms of an isolated single adatom shows a singularity at E = 0 [Fig. 3(c)]. In contrast, an enhancement of the LDOS in a broad energy region -t &lt; E &lt; t is observed on certain carbon atoms in the vicinity of the adatom dimer [Fig. 3(d)]. The numerical LDOS for selected impurities compared to the results of analytical Green's function calculations shows no discrepancies [36].

Based on these observations we conclude that the residual peak at E = 0 in the case of the correlated distribution is due to the n &gt; 2 adatom aggregates, which break the bipartite symmetry of graphene (that is, populate different number of sites in the two sublattices). This gives rise to resonant modes at E = 0 akin to isolated adatoms. All oddn aggregates and certain configurations of evenn clusters lead to resonant modes at E = 0 . Judging on the cluster size distributions [Fig. 2(d)], the largest contribution comes from adatom trimers. This suggests that adatom aggregation has strong effects on the electronic transport properties that we investigate by performing a scaling analysis of conductivity g using the Landauer-Büttiker approach [37]. In this approach, the conductance G E ( ) is given as G E ( ) = G T E 0 ( ) , where G 0 is the conductance quantum and T E ( ) is the transmission probability across the scattering region at energy E . We assume a two-terminal device configuration with a scattering region of width W = 40 nm perpendicular to the current direction, and of variable length 1 nm &lt; L &lt; 60 nm. The scattering region is attached to pristine graphene contacts and populated by adatoms according to concentrations and statistical distributions discussed above. Further details of our methodology are given in the Supplementary Material [31].

The characteristic functional laws for the conductivity g = G × L/W in the ballistic, diffusive and localized transport regimes are g ∝ L , g = const and g ∝ exp( -L/ξ loc ) , respectively, where ξ loc is the localization length. In the localized regime ln g ( ) follows a broad positively-skewed distribution, which means that g can show strong fluctuations depending on the exact configuration of defects, especially in the presence of strong localization [38]. An estimate of the mean value for such a distribution is given by the typical conductivity g typ = exp ln( 〈 g ) 〉 [38, 39]. In our scaling analysis, g typ has been obtained averaging over 9600 disorder realizations. Figure 4(a,b) shows g typ as a function of scattering region length L at different energies for the random and correlated impurity distributions, both at x = 5% concentration. We observe a short transition from ballistic to diffusive and subsequently to localized regime within the first 10 nm. The crossover lengths are expected to be of the order of the elastic mean free path ξ el and localization length ξ loc , respectively. The general trend is that the localized regime is accentuated at low energy, whereas at higher energy the onset of exponential decay occurs at larger L and the absolute slope of the conductance curves is smaller. The scaling of g also depends strongly on the impurity concentration. For low adatom concentrations ( x &lt; ∼ 0 5% . ) the onset of the localized regime is only observable in the vicinity of the Dirac energy (see Supplementary Material [31] for complete results).

We stress that the conductance curves vary smoothly in the whole range of E and x , never showing singularities which would indicate a phase transition such as the metal-insulator transition (MIT). Thus, we ascribe the non-observance of the localized regime to an insufficient scattering region length of our model, which is shorter than ξ loc for some choices of E and x . From Fig-

Figure 4: (Color online) (a,b) Averaged conductivity g typ for the random and correlated adatoms distributions, respectively, as a function of scattering region length L at different energies E and at x = 5% concentration. (c,d) Conductivity g as a function of energy E calculated using the kernel polynomial method for both adatom distributions at concentration x = 0 5% . and x = 5 0% . , respectively. The dashed lines show the minimum conductivity g = 4 /π .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0245v1-figure-ee25e3f8-97b4-4e3c-9f6f-4436ea372e8a.jpg)

ure 4(a,b) it follows that the presence of spatial correlation between adatoms enhances the conductance by up to five orders of magnitude in the vicinity of the Dirac point ( E = 2 7 . meV). This is a direct consequence of the suppressed weight of low-energy resonant states, as explained above. A closely related effect is a significant increase of the localization length ξ loc at all energies upon adatom aggregation. ξ loc was obtained by fitting conductance curves to the expected law g typ ∝ exp( -L/ξ loc ) [40], and is shown in Fig. 5(a). At x = 5% , the localization length is well defined for the entire range of investigated energies -1 eV &lt; E &lt; 1 eV, whereas at lower concentrations it is well defined only in proximity of the Dirac point (see Supplementary Material [31]). However, as long as a finite positive ξ loc can be determined, it proves to be up to an order of magnitude larger for the correlated adatom distribution compared to the random case.

Our Landauer-Büttiker results are complemented by the Kubo formula calculations based on the kernel polynomial method (KPM) [41]. This approach allows for a direct calculation of the DOS and conductivity independent of sample geometry and contacts. Figure 4(c,d) shows the conductivity g for adatom concentrations x = 0 5% . and x = 5% , and both random and correlated adatom distributions. In comparison to the LandauerBüttiker ansatz, the results are similar at both concentrations, confirming the enhancement of conductivity upon aggregation of adatoms. The only difference appears when the charge-carrier energy E falls into the resonant

Figure 5: (Color online) (a) Localization length ξ loc as a function of charge-carrier energy E for the case of random and correlated adatom distribution at x = 5% . (b) Charge-carrier localization length ξ loc at low energy ( E = 10 -3 t = 2 7 . meV) as a function of simple concentration x and effective concentration x glyph[star] for the random and correlated distributions.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0245v1-figure-90f6b7f0-e677-48d1-aba2-1c0eee713041.jpg)

peak region where the KPM conductivity exceeds the Landauer-Büttiker results. This can be proved to be an effect of the increased DOS due to the formation of an impurity band. While the KPM conductivity reflects the excess of states at the Fermi energy, in the LandauerBüttiker ansatz the number of charge carriers is limited by the DOS of the pristine graphene leads, and subsequently the conductivity is determined only by the mean free path and the localization length of the scattering region. We note that upon doping the leads, the LB conductivity at low energies increases while the localization length is unaffected (see Supplementary Material [31]).

On the basis of the identification of resonant adatom clusters we introduce an effective concentration

$$x ^ { * } = \frac { 1 } { N _ { C } } \sum _ { i } N _ { i } | n _ { i } ^ { A } - n _ { i } ^ { B } |,$$

where N i is the number of instances of the cluster configuration i , and n A i ( n B i ) is the number adatoms bound to carbon atoms in sublattice A(B) in this configuration. The main contribution to x glyph[star] comes primarily from adatom trimers as discussed above. We assume that the contribution of non-resonant clusters to the total scattering cross-section can be neglected. It should be stressed that this approach accounts for the diverse scattering effect of different clusters. In order to validate the applicability of effective concentration x glyph[star] , we compare ξ loc ( x ) with ξ loc ( x glyph[star] ) at low charge-carrier energy for the correlated impurity case [Fig. 5(b)]. One can see that replacing x with effective concentration x glyph[star] brings ξ loc ( x glyph[star] ) evaluated for the case of the correlated impurity distribution in good agreement with ξ loc ( x ) calculated for the random distribution. The agreement is particularly good in the low concentration regime where randomly distributed impurities consist almost exclusively of isolated monomers. Deviations at higher adatom concentrations can be ascribed to the inter-cluster interference effects, which become important at the reduced average distance between

the clusters.

To summarize, hydrogen adsorbed on graphene has a strong tendency toward aggregation, resulting in the formation of small clusters and fully eliminating isolated adatoms. Some of the larger clusters, notably trimers, are responsible for the residual resonant scattering, but the overall conductance and localization length dramatically increases upon aggregation. Within the range of parameters investigated in our work, we find no metal-insulator transition, with the graphene spectrum being fully localized. The predicted effects of adatom aggregation can be investigated experimentally by varying the temperature regimes, since the diffusion of hydrogen adatoms occurring at normal conditions can be effectively suppressed at low temperatures. Alternatively, time-resolved transport measurements should evince a rise in conductivity upon formation of adatom clusters.

We would like to thank C. W. J. Beenakker, P. Lugan and V. Savona for discussions. This work was supported by the Swiss National Science Foundation (Grant No. PP00P2\_133552).

- [1] A. K. Geim and K. S. Novoselov, Nature Mater. 6 , 183 (2007).
- [2] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, Rev. Mod. Phys. 81 , 109 (2009).
- [3] S. Das Sarma, S. Adam, E. H. Hwang, and E. Rossi, Rev. Mod. Phys. 83 , 407 (2011).
- [4] J. Tworzydło, B. Trauzettel, M. Titov, A. Rycerz, and C. W. J. Beenakker, Phys. Rev. Lett. 96 , 246802 (2006).
- [5] M. I. Katsnelson, K. S. Novoselov, and A. K. Geim, Nature Phys. 2 , 620 (2006).
- [6] Y. Zhang, Y.-W. Tan, H. L. Stormer, and P. Kim, Nature (London) 438 , 201 (2005).
- [7] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, M. I. Katsnelson, I. V. Grigorieva, S. V. Dubonos, and A. A. Firsov, Nature (London) 438 , 197 (2005).
- [8] T. O. Wehling, S. Yuan, A. I. Lichtenstein, A. K. Geim, and M. I. Katsnelson, Phys. Rev. Lett. 105 , 056802 (2010).
- [9] F. Gargiulo and O. V. Yazyev, Nano Lett. 14 , 250 (2014).
- [10] Z. H. Ni, L. A. Ponomarenko, R. R. Nair, R. Yang, S. Anissimova, I. V. Grigorieva, F. Schedin, P. Blake, Z. X. Shen, E. H. Hill, K. S. Novoselov, and A. K. Geim, Nano Lett. 10 , 3868 (2010).
- [11] T. O. Wehling, A. V. Balatsky, M. I. Katsnelson, A. I. Lichtenstein, K. Scharnberg, and R. Wiesendanger, Phys. Rev. B 75 , 125425 (2007).
- [12] B. Kramer and A. MacKinnon, Rep. Prog. Phys. 56 , 1469 (1993).
- [13] F. Evers and A. D. Mirlin, Rev. Mod. Phys. 80 , 1355 (2008).
- [14] A. Bostwick, J. L. McChesney, K. V. Emtsev, T. Seyller, K. Horn, S. D. Kevan, and E. Rotenberg, Phys. Rev. Lett. 103 , 056404 (2009).
- [15] Y. Song, H. Song, and S. Feng, J. Phys.: Condens. Mat-

ter 23 , 205501 (2011).

- [16] R. Jayasingha, A. Sherehiy, S.-Y. Wu, and G. U. Sumanasekera, Nano Lett. 13 , 5098 (2013).
- [17] S. Adam, S. Cho, M. S. Fuhrer, and S. Das Sarma, Phys. Rev. Lett. 101 , 046404 (2008).
- [18] N. Shon and T. Ando, J. Phys. Soc. Jpn. 67 , 2421 (1998).
- [19] T. O. Wehling, M. I. Katsnelson, and A. I. Lichtenstein, Phys. Rev. B 80 , 085428 (2009).
- [20] A. Ferreira, J. Viana-Gomes, J. Nilsson, E. R. Mucciolo, N. M. R. Peres, and A. H. Castro Neto, Phys. Rev. B 83 , 165402 (2011).
- [21] A. Cresti, F. Ortmann, T. Louvet, D. Van Tuan, and S. Roche, Phys. Rev. Lett. 110 , 196601 (2013).
- [22] J. J. Palacios, J. Fernández-Rossier, and L. Brey, Phys. Rev. B 77 , 195428 (2008).
- [23] N. Leconte, D. Soriano, S. Roche, P. Ordejon, J.-C. Charlier, and J. J. Palacios, ACS Nano 5 , 3987 (2011).
- [24] G. Trambly de Laissardière and D. Mayou, Phys. Rev. Lett. 111 , 146601 (2013).
- [25] D. W. Boukhvalov, M. I. Katsnelson, and A. I. Lichtenstein, Phys. Rev. B 77 , 035427 (2008).
- [26] Y. Lin, F. Ding, and B. I. Yakobson, Phys. Rev. B 78 , 041402 (2008).
- [27] O. V. Yazyev, Phys. Rev. Lett. 101 , 037203 (2008).
- [28] C. P. Herrero and R. Ramírez, Phys. Rev. B 79 , 115429 (2009).
- [29] M. Moaied, J. A. Moreno, M. J. Caturla, and J. J. Palacios, arXiv:1405.3165 (2014).
- [30] O. V. Yazyev and L. Helm, Phys. Rev. B 75 , 125408 (2007).
- [31] See Supplemental Material at [URL will be inserted by publisher] for detailed description of the methodology and additional calculations.
- [32] In the case of hydrogen adsorption on both sides of graphene the fitted parameters γ 1 = -1 461 . eV and γ 2 = 0 342 . eV. The resulting cluster size distributions, however, are very close to those obtained for single-side adsorption [31].
- [33] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, J. Chem Phys 21 , 1087 (1953).
- [34] V. M. Pereira, F. Guinea, J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H. Castro Neto, Phys. Rev. Lett. 96 , 036801 (2006).
- [35] D. M. Basko, Phys. Rev. B 78 , 115432 (2008).
- [36] F. Gargiulo, G. Autès, and O. V. Yazyev, to be published (2014).
- [37] M. Büttiker, Y. Imry, R. Landauer, and S. Pinhas, Phys. Rev. B 31 , 6207 (1985).
- [38] D.-H. Choe and K. J. Chang, Nano Lett. 12 , 5175 (2012).
- [39] A. Uppstu, Z. Fan, and A. Harju, Phys. Rev. B 89 , 075420 (2014).
- [40] P. A. Lee and T. V. Ramakrishnan, Rev. Mod. Phys. 57 , 287 (1985).
- [41] A. Weiße, G. Wellein, A. Alvermann, and H. Fehske, Rev. Mod. Phys. 78 , 275 (2006)

## Supplementary Material for 'Electronic Transport in Graphene with Aggregated Hydrogen Adatoms'

## A. First-principles calculations of the energies of hydrogen adatom clusters

First-principles calculations of the interaction energies of hydrogen adatoms on graphene have been performed within the density functional theory (DFT) framework employing the generalized gradient approximation (GGA) to the exchange-correlation functional [S1]. Ultrasoft pseudopotentials [S2] for carbon and hydrogen atoms have been used in combination with a plane-wave basis set with a kinetic energy cutoff of 30 Ry for the wavefunctions. Models of hydrogen adatom clusters are based on a graphene 6 × 6 supercell with 15 Å of vacuum separating the periodic replicas. We used a 2 × 2 × 1 Monkhorst-Pack k-point mesh for the Brillouin zone integration [S3]. All hydrogen adatom cluster models were relaxed until a maximum force of 0.15 eV/Å on individual atoms was reached. We verified that the chosen parameters provide sufficiently accurate total energies. All calculations have been performed using the pwscf code of the Quantum ESPRESSO package [S4].

The interaction energy of a cluster of hydrogen adatoms calculated from first principles E DFT is defined as

$$E _ { \text{DFT} } = E _ { \text{gr+nH} } - E _ { \text{gr} } - n \left ( E _ { \text{gr+H} } - E _ { \text{gr} } \right ),$$

where E gr+ n H , E gr and E gr+H are the total energies of graphene with a cluster of n hydrogen adatoms, pristine graphene and graphene with a single hydrogen adatom, respectively.

In the main text of our manuscript we focused on the situation where hydrogen adatoms are deposited on a single side of graphene. The values of the fitted parameters in expr. (1) of the main text are γ 1 = -1 182 . eV and γ 2 = 0 484 . eV. This scenario is relevant to the case of graphene on a substrate, however in the situation of suspended graphene both sides of graphene are available for binding adatoms. We investigated this situation by studying the same set of adatom clusters as shown in Fig. 1(b) of the main text, but with adatoms placed on the opposite sides of graphene sheet when functionalized carbon atoms belong to different sublattices. The fitted interaction parameters are γ 1 = -1 461 . eV and γ 2 = 0 342 . eV. The excellent agreement between the estimated interaction energies ˜ E and the first-principles values E DFT is illustrated in Fig. S1(a). This case is characterized by a stronger attractive contribution and a weaker repulsion compared to the single-side adsorption, thus reflecting the known tendency of forming more stable adatom aggregates upon adsorption on both sides [S5, S6]. The cluster size distributions P n ( ) calculated for single-side and both-sides adsorption at adatom concentration x = 5% and T = 300 K are compared in Fig. S1(b). While the distributions are qualitatively very similar, one notes that both-sides adsorption exhibits a somewhat stronger tendency to form larger clusters.

(a)

n

Figure S1: (Color online) (a) Predicted energy ˜ E of aggregation of hydrogen adatoms as a function of aggregation energy E DFT calculated from first principles for the set of small clusters shown in Fig. 1(b) of the main text with adatoms adsorbed on both sides of graphene. (b) Comparison of the cluster size distributions P n ( ) for the cases of single-side and both-sides adatom adsorption at x = 5% concentration and T = 300 K.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0245v1-figure-547b841d-fc1a-4b2b-8329-87802102306c.jpg)

## B. Monte-Carlo simulations of the hydrogen adatoms aggregation

We performed Monte-Carlo simulations with an elementary trial move being the displacement of a randomly chosen adatom to a random carbon atom not populated by another adatom. This move insures the fulfillment of detailed balance. The Metropolis algorithm has been employed for the acceptance/rejection criterion. Once a move has been performed, the system is updated from the old configuration S old to the new one S new with a probability

$$P ( S _ { \text{old} } \rightarrow S _ { \text{new} } ) & = \min ( 1, e ^ { - \beta [ \tilde { E } ( S _ { \text{new} } ) - \tilde { E } ( S _ { \text{old} } ) ] } ),$$

where β is the inverse temperature 1 / k ( B T ) . In all our simulations T = 300 K . The number of equilibration steps N eq disregarded from statistical sampling varied between 10 6 and 1 6 . × 10 7 , depending on adatom concentration (larger concentrations need more equilibration steps). The total number of steps in our simulations varied between 10 7 and 3 × 10 7 . The equilibration efficacy of our simulation was tested by comparing certain equilibrium properties such as the total energy and cluster size distribution. The properties obtained from Monte-Carlo simulations performed with and without temperature annealing show negligible differences.

## C. Landauer-Büttiker electronic transport calculations and scaling analysis of conductivity

In order to investigate the transport properties of graphene with resonant scattering impurities we perform Landauer-Büttiker calculations in a two-terminal configuration with a scattering region composed of hydrogenated graphene attached to two semi-infinite leads of pristine graphene, as shown in Fig. S2. The overall configuration is periodic along the transverse direction y . For such a setup, conductance as a function of energy G E ( ) is given by

$$G ( E ) = \frac { W } { 2 \pi } \int _ { - \tilde { w } } ^ { \tilde { w } } T ( E, k _ { \| } ) d k _ { \| },$$

Figure S2: (Color online) Schematic drawing of the two-terminal configuration employed for investigating the transport properties of graphene with hydrogen adatoms. The transport direction is along the x axis while the system is periodic along the y axis. The periodicity along the y axis is W = 40 nm. The unit cells of the left and right leads composed of pristine graphene are indicated by blue lines. The scattering region is populated by adatoms either randomly or according to the configurations produced by Monte-Carlo simulations in the case of correlated adatom distributions. One of the principal layers (PL) in the scattering region is indicated by means of green dashed lines.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0245v1-figure-1b75fb5f-3782-4785-b4e8-b6e7c74c7438.jpg)

where T E,k ( ‖ ) is the transmission probability and k ‖ is the momentum along y [S7]. Due to the large width W = 40 nm of the model employed, transmission is only evaluated at the Γ point ( k ‖ = 0 ).

In order to calculate T E ( ) we decompose the scattering region into principal layers (PLs), that is, the layers in which the atoms are coupled at most to those located in the next layer. Since our tight-binding model is limited to first-nearest-neighbor interactions and the transport direction is oriented along zig-zag direction [see Fig. S2], the minimal width of the PL is d PL = √ 3 2 d CC , where d CC = 1 42 . Å is the carbon-carbon bond length. The Hamiltonian restricted to the i -th principal layer is H i , while t i is the tight-binding hopping matrix connecting i -th and i +1 -th principal layers. We introduce an imaginary cleavage plane between the n -th and n +1 -th principal layers dividing the system into two independent parts. We define g L n and g R n +1 as the surface Green's functions of the two non-interacting semi-infinite systems located on the left and on the right sides of this cleavage plane, respectively.

Following Ref. [S8], the transmission is given by

$$T ( E ) = \text{Tr} [ T _ { n } I m ( g _ { n } ^ { L } ) T _ { n } ^ { \dagger } I m ( g _ { n + 1 } ^ { R } ) ],$$

where operator T n is defined as

$$T _ { n } = t _ { n } ( 1 - g _ { n + 1 } ^ { R } t _ { n } ^ { \dagger } g _ { n } ^ { L } t _ { n } ) ^ { - 1 }.$$

The choice of the position of the cleavage plane is immaterial because of current conservation.

Surface Green's functions g L n and g R n can be related to the preceding (successive) surface Green's functions g L n -1 ( g R n +1 ) by applying the Dyson equations [S9]

$$g _ { n } ^ { L } = ( E - H _ { n } - t _ { n - 1 } ^ { \dagger } g _ { n - 1 } ^ { L } t _ { n - 1 } ) ^ { - 1 }$$

and

$$g _ { n } ^ { \text{R} } & = ( E - H _ { n } - t _ { n } g _ { n + 1 } ^ { \text{R} } t _ { n } ^ { \dagger } ) ^ { - 1 }.$$

Further iterations of Eqns. (S6) and (S7) reduce the problem to the knowledge of the Green's functions at the surfaces separating the scattering region from the left and right leads, g LL and g RL , that we calculated according to the analytic closed form solution described in Ref. [S10].

The time complexity of the Green's function calculation for each lead with respect to the number N lead of orbitals in the lead unit cell is O N ( 3 lead ) [S9]. On the other hand, as follows from Eqns. (S6) and (S7), the complexity of the addition of the layers required to reach the cleavage plane is O M ( × N 3 layer ) , where M and N layer are the number of principal layers and the number of orbitals in each layer, respectively. Consequently, the overall complexity of the method is cubic with respect to the width and linear with respect to the length of the system.

In order to perform our conductance scaling analysis we vary the length of the scattering region in the range L = 1 ... 60 nm by steps of 8 PLs, which corresponds to ∆ L ≈ 1 nm. At each step the right lead is moved rightwards whereas the left lead is kept fixed [Fig. S2]. The values of conductance G E ( ) are averaged over an ensemble of N ens =9600 disorder realizations for proper statistical sampling.

## D. Complete account of the results of calculations of conductivity and localization length

Conductivity g for the entire investigated range of concentrations is presented in Fig. S3(a,b). The overall enhancement of conductivity upon the aggregation of adatoms is a common feature at all investigated concentrations. This is particularly visible in the strong localization regime, that is at low energies E and large scattering region lengths L . Localization length ξ loc can only be determined for the g L ( ) curves which exhibit a well-defined negative slope in the large length region. This is the case for ξ loc &lt; 60 nm, which is the maximum scattering region length considered in our study. This does not imply that the system does not undergo localization, but rather that a longer scattering region is needed in order to estimate ξ loc correctly. For this reason many values of localization length ξ loc at x &lt; 5% , especially in the case of correlated impurities, are missing in Fig. S3(c).

## E. Effect of the doping of leads in the Landauer-Büttiker calculations

Figure S4(a,b) shows the conductivity curves g L ( ) for x = 5% concentration of randomly distributed adatoms obtained by shifting the charge neutrality point of leads by ∆ E L = -0 5 . eV and ∆ E L = -1 0 . eV, respectively.

Figure S3: (Color online) Scaling analysis of conductivity g and localization length ξ loc for concentrations x = 0 25% . ... 10% . (a,b) Conductivity g as a function of scattering region length L calculated for graphene with random and correlated adatom distributions, respectively, at charge-carrier energies 0 eV &lt; E &lt; 1 eV. (c) Localization length ξ loc as a function of charge-carrier energy E for random and correlated adatom distributions.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0245v1-figure-c8bce11e-e11d-42e9-913f-d3f4059f5992.jpg)

The effect of the doping of leads is two-fold. Firstly, the DOS of pristine graphene increases away from the charge neutrality point. Hence, upon doping the number of transport channels increases, which may result in larger values of conductance g . This is particularly important when the scattering region has an enhanced DOS at zero energy, such as graphene with resonant impurities. A comparison of Fig. S4(a) and Fig. S3(a) for x = 5% shows that conductivity is indeed enhanced at low energies with a crossing of the g L ( ) curves at L ≈ 10 nm. Secondly, doping results in a mismatch between the Fermi wavelength of the leads and that of the scattering region, which has a detrimental effect on conductance g . This effect is expected to be more pronounced at higher doping and high energy, where localization plays a smaller role. Indeed, at higher doping ( ∆ E L = -1 0 . eV) and high energy ( E = 1 0 . eV) the conductivity is reduced, notably at short distance L &lt; 10 nm [Fig. S4(a,b)]. On the other hand, at large distances a general increase of the conductance is progressively restored since the conductivity becomes predominantly determined by the localization of the wavefunction and the increased number of available states. Finally, as shown in Fig. S4(c), localization length ξ loc is practically unaffected by the doping of the leads, since it is an intrinsic property of the scattering region.

Figure S4: (Color online) (a,b) Conductivity g as a function of scattering region length L for graphene with randomly distributed adatoms at x = 5% . The charge neutrality point of the leads has been shifted by (a) ∆ E L = -0 5 . eV and (b) ∆ E L = -1 0 . eV, respectively. (c) Localization length ξ loc as a function of charge-carrier energy E for the two investigated lead doping levels compared to undoped leads ∆ E L = 0 eV.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0245v1-figure-7ac379e7-08ef-45c2-a904-2f60e6eef58a.jpg)

## F. Electronic transport calculations within the Kubo-Greenwood formalism using the kernel polynomial method

The diagonal elements of the real part of the frequency dependent conductivity tensor in linear response theory is given by the Kubo-Greenwood formula

$$R e \{ \sigma _ { \alpha \alpha } ( \omega ) \} = \frac { \pi } { V } \int d E \frac { f ( E ) - f ( E + \hbar { \omega } ) } { \omega } \times \text{Tr} \{ \delta ( E - H ) j _ { \alpha } \delta ( E - H + \hbar { \omega } ) j _ { \alpha } \}.$$

Here, f ( E ) is the Fermi function and H the single-particle tight-binding Hamiltonian of the system under consideration, whereas the vectorial component of the current operator j α is defined later. In the thermodynamic limit for zero temperature and zero frequency this equation reduces to

$$R e \{ \sigma _ { \alpha \alpha } ( 0 ) \} = \frac { \pi \hbar } { V } \text{Tr} \{ \delta ( E _ { \text{F} } - H ) j _ { \alpha } \delta ( E _ { \text{F} } - H ) j _ { \alpha } \},$$

which defines the DC-conductivity studied in this work. The evaluation of this expression is carried out by employing the kernel polynomial method (KPM), as described in detail in Ref. [S11]. In this framework, the following matrix element density is defined

$$j ( E, E ^ { \prime } ) = \frac { 1 } { V } \sum _ { n, m } \langle n | j _ { \alpha } | m \rangle \langle m | j _ { \alpha } | n \rangle \delta ( E - \hbar { \omega } _ { n } ) \delta ( E ^ { \prime } - \hbar { \omega } _ { m } ),$$

which is then expanded up to finite order M within the two-dimensional KPM using the Jackson kernel. For the studied supercells containing up to N C ∼ 10 6 carbon atoms, we used M = 1280 . The real-part of the frequency

Figure S5: (Color online) (a,b) Conductivity g as a function of charge-carrier energy E calculated using the kernel polynomial method (KPM) for the random and correlated impurity distributions at x = 0 5% . and x = 5 0% . adatom concentrations, respectively. (c,d) Conductivity g as a function of expansion order M for the random and correlated impurity distributions at x = 0 5% . and x = 5 0% . adatom concentrations, respectively. The dashed lines indicate the minimum conductivity of graphene.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0245v1-figure-5ada8685-e904-42a5-bfdc-934a55ae0675.jpg)

dependent conductivity tensor is then given by a double integration over the matrix density for arbitrary temperature and Fermi level:

$$R e \{ \sigma _ { \alpha \alpha } ( \omega ) \} = \frac { \pi } { \omega } \int _ { - \infty } ^ { \infty } d E \int _ { - \infty } ^ { \infty } d E ^ { \prime } j ( E, E ^ { \prime } ) \times [ f ( E ) - f ( E ^ { \prime } ) ] \, \delta ( \hbar { \omega } - ( E ^ { \prime } - E ) ).$$

As a result, the DC-conductivity can be expressed in terms of the diagonal elements of j ( E,E ′ ) :

$$R e \{ \sigma _ { \alpha \alpha } ( 0 ) \} = \pi \hbar { j } ( E _ { \text{F} }, E _ { \text{F} } ).$$

The calculation of transport properties requires the definition of a current operator appropriate to the employed single-particle tight-binding model. Thus, we approximate the spatial operator r in a diagonal form:

$$r ^ { \gamma } \approx \sum _ { i } r _ { i } ^ { \gamma } c _ { i } ^ { \dagger } c _ { i }.$$

Consequently, one can derive the following current operator used in this work,

$$j _ { \gamma } = - i \frac { e _ { 0 } } { \hbar } \sum _ { i j } t _ { i j } ( r _ { j } ^ { \gamma } - r _ { i } ^ { \gamma } ) c _ { i } ^ { \dagger } c _ { j },$$

where indices i and j run over atomic positions and γ denotes the vectorial component. For details of a similar derivation see e.g. Ref. [S12].

## G. Results for the calculations of conductivity using the kernel polynomial method

Figure S5(a,b) shows conductivity g as a function of charge-carrier energy E calculated using the KPM for random and correlated impurity distributions at two different adatom concentrations x = 0 5% . and x = 5 0% . . A scaling analysis of the conductivity for three charge carrier energies ( E = 1 7 . × 10 -3 eV, E = 0 1 . eV and E = 1 0 . eV) based on the expansion order M is shown in Fig. S5(c,d) for the same values of adatom concentration x . The order M

relates to a timescale τ ( M ) [S13]. For E = 1 0 . eV and x = 5% the conductivity rises with increasing expansion order M , exhibits a maximum around M ∼ 200 and then decreases slowly again for both correlated (dashed line) and random (solid line) impurity distributions [Fig. S5(d)]. We interpret these maxima as the semiclassical values of conductivity g sc corresponding to the diffusive regime, followed by quantum corrections which result in the decrease of the conductivity. For a concentration of x = 0 5% . and E = 1 0 . eV no conductivity maximum is observed, which indicates that the results still correspond to the pre-diffusive/ballistic regime for the largest expansion order investigated. A shift of g sc to larger expansion orders for decreasing x is expected, because fewer scattering centers are present. The situation is, in general, different for small energies close to the Dirac point (here, E = 1 7 . × 10 -3 eV), in both cases of correlated and random adatom distributions as well as concentrations of x = 0 5% . and x = 5% . A distinct maximum for M &gt; 100 is completely absent, and only a reduction of the conductivity is observed for increasing M . Thus, we relate this characteristic to the quantum regime with localization effects resulting in a reduced localization length. If the expansion order is sufficiently large, the curves for E = 1 7 . × 10 -3 eV reach values well below g = 4 /π e /h ( 2 ) , and one may expect the conductivity to converge asymptotically to zero in the limit of M →∞ if the modes are completely localized.

- [S1] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77 , 3865 (1996).
- [S2] D. Vanderbilt, Phys. Rev. B 41 , 7892 (1990).
- [S3] H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13 , 5188 (1976).
- [S4] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. D. Corso, S. d. Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcovitch, J. Phys.: Condens. Matt. 21 , 395502 (2009).
- [S5] Y. Lin, F. Ding, and B. I. Yakobson, Phys. Rev. B 78 , 041402 (2008).
- [S6] D. W. Boukhvalov, M. I. Katsnelson, and A. I. Lichtenstein, Phys. Rev. B 77 , 035427 (2008).
- [S7] M. Büttiker, Y. Imry, R. Landauer, and S. Pinhas, Phys. Rev. B 31 , 6207 (1985).
- [S8] J. Mathon and A. Umerski, Phys. Rev. B 63 , 220403 (2001).
- [S9] A. Umerski, Phys. Rev. B 55 , 5266 (1997).
- [S10] S. Sanvito, C. J. Lambert, J. H. Jefferson, and A. M. Bratkovsky, Phys. Rev. B 59 , 11936 (1999).
- [S11] A. Weiße, G. Wellein, A. Alvermann, and H. Fehske, Rev. Mod. Phys. 78 , 275 (2006).
- [S12] J. M. Tomczak and S. Biermann, Phys. Rev. B 80 , 085117 (2009).
- [S13] S. Roche, Phys. Rev. B 59 , 2284 (1999)
Fernando Gargiulo Gabriel Autès Naunidh Virk Stefan Barthel Malte Rösner Lisa R. M. Toller Tim O. Wehling Oleg V. Yazyev
Electronic Transport in Graphene with Aggregated Hydrogen Adatoms
Hydrogen adatoms and other species covalently bound to graphene act as
resonant scattering centers affecting the electronic transport properties and
inducing Anderson localization. We show that attractive interactions between
adatoms on graphene and their diffusion mobility strongly modify the spatial
distribution, thus fully eliminating isolated adatoms and increasing the
population of larger size adatom aggregates. Our scaling analysis shows that
such aggregation of adatoms increases conductance by up to several orders of
magnitude and results in significant extension of the Anderson localization
length in the strong localization regime. We introduce a simple definition of
the effective adatom concentration $x^\star$ , which describes the transport
properties of both random and correlated distributions of hydrogen adatoms on
graphene across a broad range of concentrations.