Graphene Multi-Protonation: a Cooperative Mechanism for Proton Permeation
Massimiliano Bartolomei, Marta I. Hernández, José Campos-Martınez, Ramón Hernández Lamoneda
GGraphene Multi-Protonation: a Cooperative Mechanism forProton Permeation
Massimiliano Bartolomei ∗ , Marta I. Hern´andez, and Jos´e Campos-Mart´ınez Instituto de F´ısica Fundamental, Consejo Superior de InvestigacionesCient´ıficas (IFF-CSIC), Serrano 123, 28006 Madrid, Spain
Ram´on Hern´andez Lamoneda
Centro de Investigaciones Qu´ımicas, Universidad Aut´onomadel Estado de Morelos, 62210 Cuernavaca, Mor. M´exico (Dated: January 10, 2019)The interaction between protons and graphene is attracting a large interest dueto recent experiments showing that these charged species permeate through the 2Dmaterial following a low barrier ( ∼ + bond fromone to the other side of the carbon layer) and previous studies have found so farthat the energy barriers (around 3.5 eV) are too high to explain the experimentalfindings. Contrarily to the previously adopted model assuming an isolated proton,in this work we consider protonated graphene at high local coverage and explore therole played by nearby chemisorbed protons in the permeation process. By meansof density functional theory calculations exploiting large molecular prototypes forgraphene it is found that, when various protons are adsorbed on the same carbonhexagonal ring, the permeation barrier can be reduced down to 1.0 eV. The relatedmechanism is described in detail and could shed a new light on the interpretation ofthe experimental observations for proton permeation through graphene. PACS numbers: ∗ Corresponding author, e-mail:maxbart@iff.csic.es a r X i v : . [ c ond - m a t . m e s - h a ll ] J a n
1. INTRODUCTION
This work is partly motivated by recent experimental work on the permeation at roomtemperature of protons through graphene[1]. In that study, Geim and coworkers reportedconductance and mass spectroscopy measurements of proton transport through pristinegraphene and found that, in a temperature range of 270-330 K, the process exhibits anArrhenius-type behavior with a rather low activation energy (about 0.8 eV). Moreover, pro-tons permeate through this two-dimensional crystal about ten times faster than deuterons[2].These discoveries were absolutely unexpected since graphene was believed to be completelyimpermeable to all atoms and molecules under ambient conditions[3, 4]. A number of workshave subsequently appeared -both experimental[5–10] and theoretical[11–16]- not only stim-ulated by the promise of important applications in hydrogen technology but also with theaim to uncover the microscopic mechanisms underlying these observations.Despite much recent progress on the theoretical insight into the processes leading to sucha facile permeation of protons through graphene [11–16], we consider that a complete andsatisfactory understanding has not been achieved yet (e.g. see Refs.[11, 12] for some dis-cussions). The possible role of surface defects in the transport process[5, 6] has been alsoindicated. Some works[15, 16] have emphasized the role of quantum tunneling in effectivelylowering the energy barrier and on isotope selectivity, using models that assume that pro-tons/deuterons are free particles. In the experiments, however, protons are initially movingwithin an aqueous medium (hydrated Nafion or HCl solution), so other works[12–14] havemore realistically considered protons to be bound to water (as H O + ) in the reactants states.As Shi et al[12] indicate, two possible modes for the penetration of a hydrated proton canbe distinguished. The first one (dissociation-penetration) involves the removal of the protonfrom the water network as it crosses the graphene membrane through the hollow of a carbonring. It has been found that this mechanism is unlikely due to the large proton affinity ofwater[12, 13]. Interestingly, Feng et al [14] have recently found, by means of first principlescalculations, that the proton transport can be largely facilitated if various carbon atoms ofthe graphene layer (close to the crossing region) are in sp configurations due to hydrogena-tion. In the second mechanism (adsorption-penetration), a proton is transferred from theaqueous media to the graphene surface, where it becomes chemisorbed at the top of a carbonatom and, in a second step, it flips through the hollow from the original chemisorbed siteto a related one on the other side of the layer. Previous works[14] indicate that the barrierfor proton transfer from water to a chemisorbed state is small and therefore the first stepappears to be feasible at ambient conditions. However, the chemisorbed state is very stable,hence the barrier for the proton flipping becomes extremely high ( ∼ et al [21] on the electrochemical storage of hydrogen within a narrow single-walledcarbon nanotube (SWCNT), where a low energy reaction path was identified for a hydrogenatom flipping from the external to the internal side of the nanotube. In addition, it is alsoworth noting that our study apparently shows similarities with that of Feng et al [14], inthe sense that in both cases it is concluded that the barrier for proton penetration notablydecreases upon local protonation/hydrogenation. However, the involved microscopic mech-anisms are completely different. A detailed comparison with the above mentioned works isprovided at the end of Section 3.The paper is organized as follows. Computational methods are described in the followingsection. Next, results are presented and discussed, starting with an analysis of the stability ofprotonated graphene and continuing with a study of the permeation process as a function ofthe number of chemisorbed protons along a carbon ring. The report closes with a summarywhere further lines for research are indicated.
2. COMPUTATIONAL METHODS
We have carried out electronic structure calculations to study the role played by anincreasing number of protons ( n = 1 − H ) is sufficiently adequatefor the present purposes. A significant advantage of using molecular prototypes is that anarbitrary number of protons can be included in the calculation, whereas periodic calculationssuffer from the problem of having to neutralize the unit cell to converge and this procedurebecomes less reliable when increasing the number of charges.DFT calculations have been performed for the optimization of the protonated circum-coronene structures by using the PBE[26] functional together with the cc-pVTZ[27] basisset. These calculations were found to be in good agreement with benchmark MP2/aug-cc-pVTZ computations carried out for a smaller carbon plane prototype such as coronene.Additional calculations involving larger prototypes (C H and C H ) were carried outwith a reduced 6-311+G[28] basis set. We have verified that, in the case of circumcoronene,this smaller set provides energy variations that are in good agreement (within few percents)with those obtained with the cc-pVTZ[27] basis set. All reported energies correspond to sta-tionary points whose correct nature has been verified by carrying out harmonic frequencycalculations, used in turn to estimate zero-point energy and thermal corrections (at 298 Kand 1 atm) to thermodynamic properties. The enthalpy of the proton in the gas phase havebeen estimated to be RT as a result of the application of the standard statistical mechanicsand ideal gas expressions. Intrinsic reaction coordinate calculations have been employed tocheck that reactants and products are indeed connected with the transition states for variousof the permeation processes.All DFT computations have been performed by using the Gaussian 09 code[29]. n, chemisorbed protons -4-202468 p r o t on a ff i n i t y , e V -32-28-24-20-16-12-8-40 ∆ H , e V C H C H C H C H C H C H n=6n=5 FIG. 1: Sequential protonation of a graphenic single ring. Upper panel: corresponding enthalpyvariation with respect to the unprotonated graphene molecular prototype ( n =0) and isolated pro-tons. Lower panel: proton affinity of each protonated graphene prototype ((C x H y + n ) n + , n=0-5).The dashed, solid and dotted lines correspond to the C H (circumcoronene), C H (circum-circumcoronene) and C H (circumcircumcircumcoronene) prototypes, respectively.
3. RESULTS AND DISCUSSION
First, we address the question of graphene affinity for the chemical adsorption of protons.Structures of circumcoronene for the sequential sticking of protons above its central ring areillustrated in Fig. 1. In the upper panel of Fig. 1, we report the enthalpy variation (∆ H )as a function of the number n of chemisorbed protons, which corresponds to the gas phaseprocess C x H y + nH + → ( C x H y + n ) n + , (1)where C x H y represents the unprotonated ( n =0) graphene prototype. In the lower panel wedepict instead the proton affinity of each protonated graphene prototype, which is definedas -∆ H of the following process( C x H y + n ) n + + H + → ( C x H y + n +1 ) ( n +1)+ , n = 0 − . (2)It can be seen that in the case of circumcoronene (C H ) the consecutive addition ofa proton to the inner ring is an energetically favourable process: in fact, positive protonaffinities are obtained for n up to 3, that is up to four protons could be adsorbed. A similarresult was previously theoretically described[30] for the multiprotonation of benzene, forwhich the addition of up to 3 protons was found to lead to stable molecular structures, thatis preserving the hexagonal ring. In this case, it seems also clear that the stability of themost protonated ( n =5,6) fragments depends on the size of the considered graphene flakes.In fact, when larger graphene molecular prototypes (C H and C H ) are taken intoaccount we observe (upper panel of Fig. 1) progressively larger enthalpy variations whichseem to tend towards a converged profile as a function of n . In particular, in the case of thecircumcircumcircumcoronene (C H ) prototypes, positive proton affinities are found forall considered protonated fragments, that is the chemisorption of up to 6 protons is feasibleand leads to stable molecular structures. These results suggest that the saturation of agraphenic ring with protons is energetically possible in the gas phase and we can expecteven more favorable proton affinities for larger graphene flakes.It is also crucial for the adopted model to mention that the local properties of the multi-protonated site are found to barely depend on the size of the graphene flake. For example,partial charges and bond distances of the central ring of a 6 times protonated circumcoroneneare very similar to those of an analogously protonated circumcircumcoronene. Note that thenet charge of this central ring is about +0.5 so that an excess charge of about +5.5 spreadsover the rest of the flake. It is reasonable that this excess charge becomes more easilydistributed as the size of the prototype increases, hence this feature must be at the originof the larger proton affinities of the bigger prototypes (Fig. 1). In addition, we have alsonoticed that, once the C-H + bond is formed, most of the proton character is lost as thepartial charge on hydrogen is very close (slightly larger) to that typical of an usual C-Hbond. Nevertheless, to remind the reader that the whole graphene flakes have the chargeof the added protons, the positive charge on the hydrogen atoms will be retained in thenotations used below.Having established the stability of the graphene prototypes when several protons arechemisorbed on a benzenic ring, we start studying the permeation process for one and twochemisorbed protons as well as the underlying microscopic mechanism. In Fig. 2, the protonpenetration from one side to the other of the carbon plane is considered for the case of a single ReactantsStateProducts(1 protonTransition(TS) n=1 n=2 through hollow TS insertion TS flipped) C a C b FIG. 2: Most favorable proton flipping mechanisms for the consecutive protonation (up to twoprotons) of a graphenic single ring. The first row shows the reactants configuration with adjacentprotons chemisorbed on the same side of the carbon plane. The second row shows the transitionstate (TS) configuration. For n =2, carbon atoms bound to the flipping proton and to the adjacentproton are indicated as C a and C b , respectively. The last row shows the products configurationwith one proton flipped on the other side of the carbon plane. The corresponding electronic energyand enthalpy balances are reported in Table I. chemisorbed proton ( n = 1) as well as for that of two protons attached to two consecutivecarbon atoms ( n = 2). For the isolated proton the transition state (TS) corresponds to aplanar structure: the C-H + bond rotates to locate the proton near the center of the ring(termed as “hollow” TS) and continues the rotation to end up in an equivalent chemisorptionstate at the other side. Notice that at the TS the original C-H + bond has been broken.However, as pointed out by Miao and coworkers[4], there are chemical interactions of theproton with the surrounding carbon atoms in the ring. Still, the corresponding activationenergy (∆ E a ) is quite high and close to 3.5 eV as reported in Table I. Similar results have ReactantsStateProducts(1 protonTransition(TS) n=5 n=6 insertion TS flipped)
FIG. 3: As in Fig. 2 but with five and six protons chemisorbed on the graphenic single ring. already been obtained from previous calculations at the DFT level[4, 11, 14]. However, fortwo nearby chemisorbed protons the most favourable reaction path is quite different: therelated TS is not planar and we observe the insertion of the proton through the middle ofthe C-C bond connecting to the closest chemisorbed additional proton (“insertion” TS). Forthis to occur, the length of that C-C bond changes from 1.55 (reactants) to 2.39 ˚A (TS), sothe bond effectively breaks. Interestingly, this breaking of the C-C bond does not provokea large energetic penalty as could be expected. On the contrary, as shown in Table I, thecorresponding activation energy is about 2.8 eV, lower than that of the single chemisorbedproton. In part this is due to the fact that, during the insertion, the C-H + bond is preservedas opposed to the n = 1 case where such a bond is broken.To get insight into this new mechanism, we have examined in detail the structures of thecorresponding stationary points (right panel of Fig. 2). Let us denote H + a and H + b as theflipping and adjacent protons, respectively, and C a and C b , as the carbon atoms linked tothose protons (see right central panel of Fig. 2). Initially, these carbon atoms exhibit fourbonds in a sp -like hybridization (i.e., C a is bound to H + a , to C b and to two other adjacent TABLE I: Electronic energy and enthalphy variations associated to the most favorable singleproton flipping process for an increasing number ( n ) of adjacent chemisorbed protons on a grapheniccarbon ring. The cases of n =1,2,5 and 6 are illustrated in Figs. 2 and 3. The first two columnsshow the energy balances between reactants (R) and transition state (TS) while the last columns,those between reactants and products (P). The last line reports the values for the arrangementof five protons on a ring plus one proton on a surrounding ring (Fig. 5), leading to the lowestactivation energy. Values are in eV. R → TS R → Pchemisorbed ∆ E a ∆ H a ∆ E r ∆ H r protonsn=1 3.44 3.24 0.00 0.00n=2 2.77 2.65 -0.39 -0.39n=3 2.29 2.17 -0.38 -0.42n=4 1.76 1.69 -0.37 -0.37n=5 1.53 1.44 -0.46 -0.47n=6 1.61 1.50 -1.21 -1.20n=5+1 1.01 0.95 -0.55 -0.56 carbon atoms of the graphenic network, the corresponding C-C distances, 1.50-1.55 ˚A, beingtypical of single bonds). At the TS, the C a -C b bond is broken but the C a atom keeps itsbonding with the three remaining atoms, showing, however, quite different bond distancesand angles (C-C distances are reduced to 1.36 ˚A). Indeed, these four atoms nearly exhibita planar geometry (with (cid:100) C − C a − H + a and (cid:100) C − C a − C angles of 113 o and 134 o , respectively),indicating that the C a atom approximately adopts a sp configuration. Analogously, thethree remaining bonds around the C b atom have transformed towards a trigonal planargeometry (angles (cid:100) C − C b − H + b and (cid:100) C − C b − C being 118.5 o and 122 o , respectively), thusshowing also C b a sp configuration. Although the distance between the flipping proton andthe opposite carbon atom at the TS, H + a -C b , is short (1.32 ˚A) it does not correspond to abonding interaction but rather to a repulsive one. This is to be expected considering that,as stated before, the proton character in the C a H + a bond has been mostly lost and thereforeit can only experience a steric repulsion with the C b atom. This will be further confirmed0below as a correlation between the H + a -C b distances and the barrier heights. Given that a C-Cbond has been broken it is natural to enquire whether an open-shell structure with unpairedelectrons characterizes the TS. We have performed stability tests of the closed-shell PBEdeterminant which prove it is stable and corresponding to the lowest energy solution (fromthe comparison with triplet states and unrestricted solutions). Therefore, as a qualitativeexplanation for the relatively low activation energy of this insertion process, it can be saidthat the cost of breaking the C-C bond at the TS is compensated by the transformationof the remaining bonds to stronger sp -type ones. We have checked that, when additionaladsorbed protons are added to the same ring of carbon atoms, the penetration process mostfavourably occurs through the same kind of insertion TS, as can be seen, for instance, inFig. 3 for the case of five and six protons.Additional adsorbed protons have been consecutively added along the central ring of thecircumcoronene molecule and the results for the proton permeation are collected in TableI and in Fig. 4. Note that the flipping proton is that located at one end of the row ofchemisorbed protons, as indicated schematically in the upper part of Fig. 4. It is found thatthe activation energy further decreases with the protonation degree, reaching a minimum of∆ E a =1.53 eV for n = 5 (less than half of the value for n =1) and then slightly increasingto 1.61 eV for n = 6. The corresponding enthalpies are roughly 0.10 eV smaller than theactivation energies, in accordance with the finding that the zero point energy is larger forreactants than for the TS. On the other hand, as seen in Fig. 4 (middle panel), the C a -C b distance in the reactants state monotonously increases with n up to n = 5. This is aconsequence of the amplification of the ring area, in turn due to the increasing number ofsingle C-C bonds and effects of steric repulsion between hydrogens. This feature probablyfacilitates a larger C a -C b distance in the TS (as indeed seen in the lower panel of Fig. 4) andtherefore an easier insertion of the proton between the two carbon atoms (a lower energybarrier). In fact, the C b -H + a distance increases from 1.32 ˚A for n =2 to 1.52 ˚A for n =6whereas the C a -H + a distance stays in the range of 1.07-1.08 ˚A in the whole n range studied.¿From Table I it can also be seen that (except for n =1 where initial and final states areequivalent) the flipping process is exothermic with electronic energy and enthalpy variations(∆ E r and ∆ H r ) ranging between -0.4 to -0.5 eV for n =2-5. This result is not surprisingsince in related systems as hydrogenated graphene[31] the most stable states involve adjacentcarbon atoms, each one linking hydrogen atoms at opposite sides of the layer. Moreover,1 ReactantsTransition State C a - C b d i s t an c e ( Å ) ∆ E a ( e V ) b+a+ Flipping proton (H ) Adjacent proton (H ) − FIG. 4: Upper panel: Activation energy ∆ E a (in eV) for the permeation of a proton through agraphene prototype (circumcoronene) as a function of the number of chemisorbed protons alongthe central carbon ring. Initial configurations of the protons are schematically depicted in theupper part of the figure, where red and dark blue circles represent flipping (H + a ) and adjacent(H + b ) protons. Middle panel: distances between the corresponding carbon atoms C a and C b in thereactant state. Lower panel: same as middle panel for the transition state. note that ∆ E r drops to about -1.20 eV for n = 6, indicating that the permeation processis globally more favorable when the graphenic ring protonation is complete. In this case asp configuration of the carbon atom bound to the flipped proton appears to be even morefeasible as it is connected with two carbon atoms supporting protons on the other side, inan arrangement that is more similar to the typical chair conformation of graphane[31].Interestingly, we have found that the permeation process exhibits different activationbarriers and exothermicities depending on the position of the proton within a given rowalong a carbon ring. For instance, we have studied the case where the flipping proton isthe central one in a row of five protons and compare the results with those of the n = 52 ReactantsStateProducts(1 protonflipped−in)Transition(TS) insertion TS n=5+1
FIG. 5: Same as in Figs.2 and 3 for the case of five protons on a ring plus another one on a nearestneighbor ring. case reported above, which shares the same reactant state but where the flipping proton isone at the end of the row (left panel of Fig. 3). The process occurs through the same kindof insertion TS but the activation energy is considerably larger (2.54 eV vs. 1.53 eV). Wehave noticed that at the TS the C a -C b distance is smaller than in the previous case; thereare also some other differences in size, shape and flatness of the central and adjacent ringsbetween the two TS being compared. In addition, exothermicity of this process is larger(-0.83 eV) than for the flipping of the external proton (-0.46 eV). As already discussed for n =6, the greater stability of the products can be related with a larger facility of the carbonatom linking the flipped proton to be arranged in a sp geometry.As can be foreseen from the previous results, it is a real challenge to identify the opti-mum configuration for the proton transport among a large number of different initial protondistributions. We have analyzed some more candidates and in Table I we show the ob-tained optimum configuration, named 5 + 1, corresponding to five protons adsorbed on thecentral ring plus another one on a nearest neighbor ring, as reported in Fig. 5. The relatedactivation energy is only 1.0 eV which already can be considered consistent with the ex-perimental determination[1], especially taking into account that a further reduction when3properly including tunneling effects could be expected. This brings strong support to thenewly proposed mechanism which could provide a new perspective for the explanation ofthe experimental findings.Moreover, we would also like to address the question of the reliability of present findingswith respect to the size of the finite graphene prototype. To do that we have calculatedthe corresponding activation energy for the proton flipping process also for circumcircum-coronene (C H ) whose estimations. Even if the energy barriers are globally larger for thelarger prototype a similar trend can be observed in both cases: in fact, they decrease withthe number of chemisorbed protons up to n =5 and, more importantly, reach a similar valuearound 1 eV for the n =5+1 case. These results suggest that low activation energies, whichare compatible with the experimental findings, are expected also for larger graphene flakes.Finally, there have been two previous theoretical studies which deserve special mentionin connection with the present work, both of them already briefly mentioned in the Intro-duction. First we recall the work of Feng et al [14] where the importance of locally saturatingcarbon chemisorption sites via hydrogenation together with the associated change in carbonhybridization proved crucial in determining lower barriers for proton transport. We wouldlike to stress some important differences with our proposed new mechanism. First, ourmodel strictly includes the addition of protons as opposed to hydrogen atoms thus lendingour approach closer to the proton permeation experimental conditions. Second and mostimportant, the permeation process reported by Feng et al involves a completely differentmechanism, as the proton is not initially chemisorbed (as in the present case) since thenearby chemisorbed sites are already saturated by hydrogenation, this unstability in theinitial state leading to an important lowering of the barrier. Moreover, proton penetratesthrough the hole of the carbon ring, in contrast with the insertion transition state presentedhere. In this way, both mechanisms are qualitatively different, but with one not excludingthe other.The process reported here does bear a resemblance to the flip-in hydrogen insertionmechanism across a (5,5) SWCNT by Lee et al [21]. By means of periodic density-functionaltight-binding calculations, they found a low energy path (with ∆ E a = 1.51 eV) for hydrogenatom insertion into the middle of a very stretched C-C bond (detailed structure of thetransition state not provided), a bond that is exothermically recovered after the hydrogenatom has flipped-in. These similarities were not necessarily foreseeable, given the differences4in the structure of the (very curved) (5,5) SWCNT and (flat) graphene and the well-knowndependence of reactivity on the curvature of the carbon substrate[32]. However, they notethat the flip-in process is efficient only if the nanotube is completely saturated, while ananalogous conclusion for protonated graphene is not suggested from the present explorations.Also, the authors report even lower activation barriers for the subsequent flipping-in ofnearby hydrogen atoms, a finding that is opposite to the test carried out with a secondproton.
4. SUMMARY AND OUTLOOK
To summarize, we have reported DFT calculations of the permeation (flipping) of a protonthrough multiprotonated graphene, using circumcoronene as its molecular prototype. A newmechanism involving relatively low energy barriers (down to about 1.0 eV) has been foundto occur when there is at least one other chemisorbed proton next to the flipping one. Thecorresponding transition state is characterized by the insertion of the proton into the middleof a C-C bond and a rearrangement of the hybridizations of the involved carbon atoms, witha recovering of that bond after the flipping is completed. The nature of this transition stateas well as the enlargement of the C-C distances in the initial state is at the origin of thereduced activation energy as compared with the flipping of an isolated chemisorbed proton.The lowest reported energy barrier is close to the experimentally measured activationenergy (0.8 eV), and some contributions not taken into account in the present work -such assolvent effects[11–13], nuclear quantum effects[13–16, 33, 34] and bias potential[1, 2]- couldactually play a role in further decreasing it. In the context of these experiments and withinthe frame of the adsorption-penetration model mentioned in the Introduction, it would beparamount to study the multiprotonation of graphene by proton transfer from the aqueousenvironment. A preliminary exploration indicates that the energy balance for the protontransfer from hydronium ions is strongly dependent on the size of the graphene prototypeand that many effects such as the formation of adducts and the simulation of the aqueousmedium should be carefully addressed. Moreover, we believe that it is worth to investi-gate whether related systems such as hydrogenated graphene or hydrogenated/protonatedhexagonal boron nitride exhibit permeation processes with mechanisms analogous to thathere reported. Work along some of these directions is in progress.5
Acknowledgments
We gratefully thank N. Halberstadt and A. Beswick for valuable directions provided. Wealso thank Maciej Gutowski for many useful discussions on extended systems. The work hasbeen funded by the Spanish grant FIS2017-84391-C2-2-P. Allocation of computing time byCESGA (Spain) is also acknowledged. RHL would like to thank CONACYT for a sabbaticalscholarship. [1] S. Hu, M. Lozada-Hidalgo, F. C. Wang, A. Mishchenko, F. Schedin, R. R. Nair, E. W. Hill,D. W. Boukhvalov, M. I. Katsnelson, R. A. W. Dryfe, et al., Nature , 227 (2014).[2] M. Lozada-Hidalgo, S. Hu, O. Marshall, A. Mishchenko, A. N. Grigorenko, R. A. W. Dryfe,B. Radha, I. V. Grigorieva, and A. K. Geim, Science , 68 (2016).[3] V. Berry, Carbon , 1 (2013).[4] M. Miao, M. B. Nardelli, Q. Wang, and Y. Liu, Physical Chemistry Chemical Physics ,16132 (2013).[5] J. L. Achtyl, R. R. Unocic, L. Xu, Y. Cai, M. Raju, W. Zhang, R. L. Sacci, I. V. Vlassiouk,P. F. Fulvio, P. Ganesh, et al., Nature Communications , 6539 (2015).[6] M. I. Walker, P. Braeuninger-Weimer, R. S. Weatherup, S. Hofmann, and U. F. Keyser,Applied Physics Letters , 213104 (2017).[7] X. Liu, Y. Gao, Y. Li, R. Wang, D. Zhu, Z. Xu, G. Ji, S. Jiang, B. Zhao, G. Yin, et al., ACSNano. , 8970 (2017).[8] M. Lozada-Hidalgo, S. Zhang, S. Hu, A. Esfandiar, I. V. Grigorieva, and A. K. Geim, NatureCommunications. , 15215 (2017).[9] S. Bukola, Y. Liang, C. Korzeniewski, J. Harris, and S. Creager, Journal of the AmericanChemical Society , 1743 (2018).[10] M. Lozada-Hidalgo, S. Zhang, S. Hu, V. G. Kravets, F. J. Rodriguez, A. Berdyugin, A. Grig-orenko, and A. K. Geim, Nature Nanotechnology. , 300 (2018).[11] J. M. H. Kroes, A. Fasolino, and M. I. Katsnelson, Phys. Chem. Chem. Phys. , 5813 (2017).[12] L. Shi, A. Xu, G. Chen, and T. Zhao, The Journal of Physical Chemistry Letters. , 4354(2017). [13] N. T. Ekanayake, J. Huang, J. Jakowski, B. G. Sumpter, and S. Garashchuk, Journal ofPhysical Chemistry C , 24335 (2017).[14] Y. Feng, J. Chen, W. Fang, E.-G. Wang, A. Michaelides, and X.-Z. Li, The Journal of PhysicalChemistry Letters. , 6009 (2017).[15] I. Poltavsky, L. Zheng, M. Mortazavi, and A. Tkatchenko, Journal of Chemical Physics ,204707 (2018).[16] J. W. Mazzuca and N. K. Haut, Journal of Chemical Physics , 224301 (2018).[17] R. Balog, B. Jorgensen, L. Nilsson, M. Andersen, E. Rienks, M. Bianchi, M. Fanetti, E. Lgs-gaard, A. Baraldi, S. Lizzit, et al., Nature Materials. pp. 315–319 (2010).[18] M. Bonfanti, S. Achilli, and R. Martinazzo, Journal of Physics. Condensed Matter (2018).[19] P. Merino, M. Svec, J. Mart´ınez, P. Jelinek, P. Lacovig, M. Dalmiglio, S. Lizzit, P. Soukiassian,J. Cernicharo, and J. Mart´ın-Gago, Nature Communications. (2014).[20] J. I. Mart´ınez, J. A. Mart´ın-Gago, J. Cernicharo, and P. L. de Andres, Journal of PhysicalChemistry C , 26882 (2014).[21] S. M. Lee, K. H. An, G. Seifert, Y. H. Lee, and T. Frauenheim, Journal of the AmericanChemical Society. , 5059 (2001).[22] M. Bartolomei, E. Carmona-Novillo, M. I. Hern´andez, J. Campos-Mart´ınez, and F. Pirani, J.Phys. Chem. C , 10512 (2013).[23] M. Bartolomei, R. P. de Tudela, K. Arteaga, T. Gonz´alez-Lezana, M. I. Hern´andez, J. Campos-Mart´ınez, P. Villarreal, J. Hern´andez-Rojas, J. Bret´on, and F. Pirani, Phys. Chem. Chem.Phys. , 26358 (2017).[24] Y. Wang, H.-J. Qian, K. Morokuma, and S. Irle, The Journal of Physical Chemistry A. ,7154 (2012).[25] E. Davidson, J. Klimes, D. Alfe, and A. Michaelides, ACS Nano. , 9905 (2014).[26] J. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. , 3865 (1996).[27] R. A. Kendall, T. H. Dunning, and R. J. Harrison, J. Chem. Phys. , 6796 (1992).[28] J. S. Binkley, J. A. Pople, and W. J. Hehre, J. Am. Chem. Soc. , 939 (1980).[29] M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman,G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson, et al., Gaussian 09 Revision E.01 ,gaussian Inc. Wallingford CT 2009.[30] R. Sumathy and E. S. Kryachko, J. Phys. Chem. A , 510 (2002). [31] M. Pumera and C. H. A. Wong, Chem. Soc. Rev. , 5987 (2013).[32] Z. Chen, W. Thiel, and A. Hirsch, ChemPhysChem , 93 (2003).[33] M. I. Hern´andez, M. Bartolomei, and J. Campos-Mart´ınez, J. Phys. Chem. A , 10743(2015).[34] A. Gij´on, J. Campos-Mart´ınez, and M. I. Hern´andez, Journal of Physical Chemistry C.121