Permeation of chemisorbed hydrogen through graphene: a flipping mechanism elucidated

Abstract The impermeability of defect-free graphene to all gases has been recently contested since experimental evidence (see Nature 579, 229–232 (2020)) of hydrogen transmission through a two-dimensional carbon layer has been obtained. By means of density functional theory computations here we elucidate a flipping mechanism which involves the insertion of a chemisorbed hydrogen atom in the middle of a C–C bond via a transition state that is relatively stable due to a sp2 rehybridization of the implicated carbon atoms. Present results suggest that transmission for hydrogenated graphene at low local coverage is highly unlikely since other outcomes such as hydrogen diffusion and desorption exhibit quite lower activation enthalpies. However, at high local coverage, with a given graphenic ring tending to saturation, the proposed flipping mechanism becomes competitive leading to a significantly exothermic process. Moreover, for a specific arrangement of four neighboring chemisorbed hydrogen atoms the flipping of one of them becomes the most likely outcome with a low activation enthalpy (about 0.8 eV), which is in the range of the experimental estimation. The effect of charge doping is investigated and it is found that electron doping can help to reduce the related activation enthalpy and to slightly enhance its exothermicity. Finally, an analysis of corresponding results for deuterium substitution is also presented.


INTRODUCTION
Two-dimensional materials are being used for developing a new class of membranes for the transport of atoms and molecules at the nanoscale [1]. As these membranes are atomically thin, permeances are expected to become much larger than for thicker materials, and nanopores can be tailored for allowing a large selectivity of the desired species with respect to other molecules in a gas or liquid mixture. Until recently, it has been believed that at ambient conditions, perfect graphene is completely impermeable to even the smallest atoms [2,3]. For example, by means of Density Functional Theory (DFT) calculations, Miao et al [4] found that protons (H + ) and atomic hydrogen (H) need to surmount rather large energy barriers (∼ 2.2 and 2.9 eV, respectively) to pass through the middle of the graphenic carbon ring. If it is assumed that these species are initially chemisorbed to graphene, even larger energy barriers were obtained for the flipping of the atoms from one side to the other of the graphene sheet (∼ 3.4 and 4.5 eV for H + and H, respectively).
Contrarily to previous predictions, in 2014 Hu et al [5] reported transport and mass spectroscopy measurements of protons permeating graphene with a much lower activation energy (∼ 0.8 eV). These findings inspired numerous studies with the aim to uncover the microscopic mechanisms underlying these observations [6][7][8][9][10][11][12][13][14], including the possible role of small defects in enhancing the proton conductivity [15][16][17]. It must be noted that, in those experiments, the graphene membrane is in contact with an aqueous medium (hydrated Nafion of HCl solution) where protons are mainly attached to water as hydronium (H 3 O + ) [8]. Removal of a proton from hydronium and a subsequent penetration of the former through graphene as a free particle seems unfeasible due to the large proton affinity to water [8,9]. However, transferring a proton from hydronium to the graphene carbon atoms leading to its chemisorption is a more likely process [10]. Recently, some of us have studied [18] a possible mechanism responsible for the flipping of protons through graphene with the assumption that they are initially chemisorbed at a high local coverage. Using DFT calculations, we found that the involved barriers notably reduce -with respect to the case of a unique chemisorbed proton (3.4 eV)-to about 1 eV when neighboring protons are adsorbed to the same carbon ring.
For two protons chemisorbed to two adjacent carbon atoms in the ring (ortho dimer), the flipping mechanism [19] involves the insertion of one proton into the middle of the C-C bond, a bond that is restored after the process is completed. Rehybridization of these carbon 4 atoms at the transition state is at the origin of the reduction of the permeation barrier, as well as the increasing saturation of the carbon ring with chemisorbed protons, a feature that enlarges the ring area for the initial state [18]. In a related study, Feng et al [10] also found that hydrogenation of graphene greatly facilitates the permeation process.
Very recently, experiments by Sun et al [20] have shown that even hydrogen, as a gas (H 2 ), permeates pristine graphene with an activation energy close to that previously found for protons. Sun et al fabricated micrometre-size monocrystalline containers tightly sealed with graphene and, upon exposure to a gas atmosphere, they measured the eventual elevation of the membrane due to permeation of the gas inside the container. In stark contrast with measurements using other gases like He, they found a noticeable permeation of H 2 through graphene, ruled by an activation energy of about 1 eV. In addition, substitution by deuterium led to an undetectable permeation of the heavier isotope, in line with the isotopic effect found for their ionic counterparts [21]. To rationalize these findings, the authors proposed that molecular hydrogen first dissociates at graphene ripples and, in a second step, the chemisorbed atoms flip to the other side but no detailed insight has been provided for the mechanism responsible for the latter step.
Strongly motivated by these novel experimental evidences, in this work we propose and elucidate a possible mechanism which accounts for the permeation of chemisorbed hydrogen atoms through graphene. The present study can also be significant for the development of new technologies for hydrogen isotopic separation [22] as well as for a deeper understanding of hydrogen intercalation, which has been evidenced in graphene/substrate interfaces [23][24][25].
Indeed, as a starting point here we consider an increasing number of hydrogen atoms chemisorbed along the rim of a given carbon ring, and transition states and energy barriers for the flipping of a hydrogen atom are determined as a function of the number of the adjacent H atoms. To this end, DFT calculations were carried out by using molecular and periodic prototypes of graphene and the competition of permeation with other relevant processes in hydrogenated graphene [32] such as diffusion, desorption and recombination are also analyzed in detail.
Finally, the effect of charge doping of the graphene support and of hydrogen substitution 5 by deuterium on the mechanism responsible for the hydrogen permeation is also studied and discussed.
The article is organized as follows. Computational methods are described in Section 2. Section 3 is devoted to the presentation and discussion of results. Finally, conclusions are given in Section 4. Some auxiliary results are reported in a Supporting Material document.

COMPUTATIONAL METHODS
DFT calculations have been performed for the optimization of the hydrogenated circumcoronene (C 54 H 18 ) and circumcircumcoronene (C 96 H 24 ) structures by using the PBE [33] functional together with the 6-311+G [34] and cc-pVTZ [35] basis sets; moreover, the density fitting method [36] has been applied to approximate the two-electron repulsion integrals. In particular, the geometry optimization have been carried out at the PBE/6-311+G level and once the stationary points have been found the corresponding energy has been evaluated at the PBE/cc-pVTZ level by performing single point calculations. The correct nature of the obtained stationary points has been verified by carrying out harmonic frequency calculations, used in turn to estimate zero-point energy (ZPE) and thermal corrections (at 298 K and 1 atm) to the estimated thermodynamic properties, namely activation and reaction enthalpies. Intrinsic reaction coordinate calculations have been employed to check that reactants and products are indeed connected with the corresponding transition states for most of the investigated cases. In order to check the reliability of the obtained results, additional calculations have been performed by using both B3LYP [37] and M062X [38] functionals to carry out estimations for the relevant thermodynamic properties for a couple of cases related to the circumcoronene finite model. These estimations are reported in Table S1 of the Supporting Material where they are compared with those carried out by using the PBE functional.
In there it can be appreciated that both absolute value and related trend of the PBE results are similar to those obtained with the B3LYP and M062X functionals, even if the latter provides in general more activated and more exothermic processes. All DFT computations have been performed by using the Gaussian 09 code [39].
In order to validate the results obtained for atomic hydrogen chemisorption by exploiting the finite model (based on C 54 H 18 and C 96 H 24 ), additional DFT calculations involving a 6 periodic model of graphene have been carried out by using the same PBE functional as implemented in the VASP code [40]: the projector augmented wave (PAW) method [41,42] has been used with a cutoff energy of 600 eV for the plane-wave basis set and the Brillouin Zone was sampled by means of a 6×8×1 Monkhorst-Pack mesh [43]. Such a large number of k-points has been chosen in order to improve the accuracy of the results: it is well known indeed that graphene electronic properties calculations require particular care in the sampling of the Brillouin zone to reach the convergence [44].
The pristine supercell sheet of graphene, constituted by 216 C atoms, has been built independently and not based on the features of the finite prototypes. To do that, a basic tetragonal (4-atoms) cell instead of the conventional hexagonal one has been used [45][46][47]. All atoms of the supercell have been fully optimized (ions and lattice parameters) till forces were < 0.01 eV/Å. The so optimized graphene layer is characterized by in-plane lattice parameters a=25.64Å and b=22.20Å. Such parameters have been kept fixed in the estimation of the reported hydrogen adsorption energies, while a sufficiently thick amount of vacuum (∼ 20Å) has been added along the non-periodic direction in order to avoid any possible spurious interaction among replicas. Also, to check the reliability of these results, we have also considered a smaller supercell (180 C atoms, optimized lateral parameters a=21.37Å and b=22.20Å) and calculated the energy variation for the single hydrogenation process, which is found to differ by only ∼0.03 eV with respect to that for the larger cell.

Chemisorption of H atoms
First, we address the question of graphene affinity for the chemical adsorption of hydrogen atoms. To this aim two different molecular graphene prototypes (circumcoronene and circumcircumcoronene) as well as a periodic model are considered and the sequential sticking of hydrogen atoms above the same carbon ring is studied. In particular, in the upper panel of Fig. 1, we report the electronic energy variation (∆E) as a function of the number n of chemisorbed hydrogen atoms (H), which corresponds to the gas phase process where C x H y represents the unsaturated (n=0) graphene prototype. In the lower panel of Fig. 1 we depict instead the electronic energy variation at each hydrogenation step, which is defined as the ∆E of the following process First of all, it can be appreciated that the obtained results do not significantly depend on shown in the case of hydrogen dimers [29], due to a nonmagnetic coupling with total spin S=0. For odd n, the two sublattices are inevitably not equally occupied leaving one electron deriving from one broken π bond unpaired: this leads to a lower energy variation for the hydrogen addition since it represents a less stable configuration with respect to the cases with even n, all characterized by paired electrons. Therefore, present results demonstrate that the saturation of a single graphenic ring with hydrogen atoms is energetically possible in the gas phase.
At this point it is useful to stress that the related hydrogen sticking barrier is low and in general theoretical estimations [32] provide values falling between 0.2 and 0.4 eV. Recently, an experimentally measurement of the hydrogen chemisorption threshold has been also obtained [48], providing a value (∼ 0.4 eV) in agreement with the theoretical estimations. Our prediction (considering the zero-point energy and thermal corrections) of the involved barrier is 0.24 eV (see Table I), which is in agreement with recent calculations obtained at a similar level of theory (See Supporting Material of Ref. [48]). The latter estimation indeed does not take into account the physisorption well [49] arising between H and the carbon support and whose inclusion should provide a better accord with the experimental chemisorption threshold. We have checked, as reported in Table I, that at each of the hydrogenation steps here considered the sticking barrier does not significantly vary for odd n, while we have been unable to find any barrier for even n. The latter is consistent with the results reported in Ref. [49] for the chemisorption of two H atoms at low temperatures: in particular in there the authors showed that for the adsorption of a second H atom the sticking barrier significantly lowers or disappears. Furthermore, the observed behavior is consistent with the closed-shell vs open-shell character explanation given above. For even n product the reaction involves the addition of two radicals, a process which is expected to be barrierless, whereas for odd n product it involves a closed-shell reactant species and the reaction has to disrupt its delocalized electron pattern, associated with a greater stability, in order to form a radical product.

Present predictions showing that a high local density of hydrogen H is both
thermodinamically and kinetically favourable seem apparently to be incompatible with the maximum experimental coverage (hydrogen/carbon ratio), which is believed to be below 40% [50,51] for atomic H and deuterium adsorbed on extended graphene. However, it must be pointed out that the detected exper- imental coverage is in general not uniform [24] over the graphene surface with both clean and high H density areas. In those areas with high coverage not only H adatom dimers [49] are found but also more complex structures due to H clustering [52], which has been shown to specially occur in curved regions and protrusions of the surface. Therefore, although the full hydrogenation of extended graphene samples has been not experimentally achieved yet as for the case of finite carbon platforms such as coronene [53], we consider that regions with a high density of H adatoms should locally exist on 2D carbon layers (see also a discussion on the topic in Section 3.3 of Ref.

Hydrogen permeation and competing processes
Having established the features of graphenic single ring being progressively saturated by H adatoms, we start studying the hydrogen permeation through the same flipping mechanism we have previously found [18] in the case of protons. In this mechanism (see upper part of This leads to a related activation energy (and enthalpy) which tends to decrease with the number of chemisorbed atoms, n, as shown in Table II and in the upper panel of Fig. 2.
As for the corresponding reaction enthalpy, it is always negative, except for n=1, and its absolute value increases with n as seen in Table II Table 1 of Ref. [18]). This therefore confirms that this flipping process is always exothermic for n >1 and becomes more favourable when the graphenic ring reaches the complete saturation. In general, it can be also appreciated that present values are slightly more exothermic (up to 0.4 eV for the highest n) than those obtained in the case of proton permeation (see Table 1 of Ref. [18]).
It must be stressed that the results in Fig. 2 obtained for circumcoronene and circumcircumcoronene are consistent and well agree to each other, confirming that both finite models for the graphene support are reliable. Therefore, for the sake of simplicity, in the following the reported results are those obtained by using the circumcoronene model.
At this is point, it is worth to consider those processes which can be competitive with that of the hydrogen insertion. Indeed, chemisorbed hydrogens can also undergo diffusion, desorption and recombination processes [32]. In fact, a sticked hydrogen could migrate to an adjacent carbon atom of the same ring; it could undergo desorption through a dehydrogenation step, the inverse process of that in Eq. 2; it could also recombine with a neighbour hydrogen atom to form one H 2 molecule, which desorbs far away. Further competitive processes, such as hydrogen migration from its chemisorption site to the nearest neighbour in the same sublattice, which is expected with a similar activation as direct diffusion, have been not taken into account for the sake of simplicity.

13
The enthalpy balances related to the considered possible outcomes for the outer chemisorbed hydrogen atom on a given graphenic ring are reported in Table III as a function of n together with sketches that schematically illustrate the involved processes. Notice that, for even n and in the case of desorption, ∆H a assumes the same value as ∆H r since for this process no transition state can be found.
The first thing to be stressed is that the hydrogen flipping through an insertion process is the only one being exothermic for all n=2-6. In fact, for both diffusion and desorption ∆H r is always positive (except n=5 for diffusion being reactant and product degenerate configurations) while recombination is the most exothermic one for n=2-5 but becomes endothermic for n=6. As for the activation enthalpy, it is clear that at all n the flipping process is that having the largest ∆H a especially for the lower n and therefore it is expected to be globally less frequent. However, for larger n the differences between the related ∆H a tend to reduce, especially for n=5. For this specific configuration the most probable event would be desorption (∆H a =0.90 eV) but the product would be less stable than the reactant and the desorbed hydrogen could stick back due to a low energy barrier (see Table I); hydrogen diffusion would almost equally occur (∆H a =1.06 eV) but it could lead to an analogous (and degenerate) configuration as that for the reactant. Also for n=5, ∆H a for the flipping process is in turn larger than those of desorption and diffusion by 0.5-0.7 eV but for the hydrogen flipped on the other side it would be less probable to flip back being the related activation enthalpy around 2.2 eV. Similar considerations could stand for the n=6, considering also that in this case ∆H r is the most negative for the flipping.
As for the recombination process, the obtained ∆H a and ∆H r for n=2 are in agreement with the theoretical estimations at a similar DFT level reported in Ref. [54]; for larger n we have reported the values related to the most favourable ("ortho" or "para" ) recombination mechanism and it can be seen that the related activation enthalphies are in general larger than those of the competing processes.
Therefore, this analysis suggests that at high coverage the hydrogen flipping through the insertion mechanism, even if not the most probable, it could occur and it would lead to the most favourable outcome (i.e. for n=6).   we have also considered alternative arrangements of a minimum number of chemisorbed atoms which could lead to a lower permeation barrier. Particularly interesting is the n=2+2 arrangement, with two nonconsecutive couples of hydrogen atoms chemisorbed on the same graphenic ring, whose possible formation path is shown in Fig. S1 of the Supporting Material together with that of the standard n=4 arrangement: in there it is shown that the n=2+2 arrangement is energetically favourable both globally and at each hydrogenation step, even if less stable than the previous n=4 case. For this arrangement it exists indeed a more elaborated mechanism which leads to the hydrogen flipping, as shown in It is worth pointing out that the C-C distance in the n=2+2 reactant configuration is about 1.62Å, which is larger than those of the previously considered arrangements and globally contributing to a less activated hydrogen insertion. In fact, in previous works devoted to the permeation of protons [10,18] it has been argued that the enlargement of the carbon ring area in the initial state (reactant) is responsible for the decrease of the activation barrier as protons are succesively chemisorbed to the ring. In the present case, an increase of strain within and around that ring is also expected due to the formation of C-C single bonds and augmented H-H steric repulsion as more atoms are chemisorbed. In particular, permeation should be facilitated if the C-C bond that is breaking at the transition state is already stretched in the reactant state. In Fig. S2 of the Supporting Material we depict the activation enthalpy as a function of the length of this C-C bond in the initial state for the various cases. As expected, this distance increases with the number of chemisorbed H atoms in the carbon ring (n=2-6) and becomes very large for the 2+2 chemisorption state. However, the activation enthalphy does not always decrease with the C-C distance, as occurred for the case of protons [18]. While the enthalpy certainly lowers along the sequence n= 3, 5 and 2+2, for n= 2, 4, and 6 it keeps high values irrespective of the elongation of the corresponding C-C bond. The latter feature must be due to the special stability of the n=even initial states, as discussed in the previous section (see also lower panel in Fig. 1). Therefore, besides the elongation of the involved bonds in the reactant state, relative destabilization of the latter (due to unpairing of electrons, for instance) is also contributing to the lowering of these permeation barriers.
In order to assess the viability of the multi-step mechanism illustrated in Fig. 3, in Table IV we also report the enthalpy balances for all possible processes originating from the n=2+2 arrangement: it can be seen that flipping is the process having in general the lowest activation enthalpy as well as a still negative reaction enthalpy. Indeed, for desorption we have reported for ∆H a the same value as for ∆H r since for this process no transition state has been obtained, as already pointed out above for the sticking barrier in the case of even n (see also Table I); therefore, if a hydrogen atom would tend to separate from the carbon plane towards desorption it could probably stick back. As for recombination, it can be seen that the formation of a H 2 molecule via an ortho mechanism is clearly more TABLE IV. Activation and reaction enthalpy for the hydrogen flipping, diffusion, desorption and recombination processes which can occur starting from the n=2+2 arrangement (see also Fig.3).
The reported ∆H a for flipping corresponds to that of the limiting step, while for recombination ∆H a and ∆H r refer to the ortho mechanism (see also

Charge doping effects
At this point it is worth noting that the corresponding permeation activation energies for protons, reported in Table 1 of Ref. [18], are in general lower than those for hydrogens (see Table II). This is particularly evident for n =4, 6 where lower values of about 1.0 eV are found in the case of protons. Therefore one may wonder whether this feature would depend on the global charge bore by the investigated system. Indeed, it has been theoretically predicted [55,56] that charge doping can actually improve the adsorption energy as well as the mobility of hydrogen atoms chemisorbed on graphene. Therefore, to investigate this effect we have considered three different charge values on the circumcoronene substrate, namely -2, 0 and +2 a.u., corresponding to electron doping, no doping and hole doping, respectively, and the related results are reported in Fig. 4. In particular, in the upper panel we report the formation enthalpy (∆H f ) of the hydrogenated carbon support with respect to the isolated support and n hydrogen atoms (see Eq. 1, with both non-hydrogenated and hydrogenated substrate bearing the charge doping). In the intermediate and lower panels of Fig. 4 we show related activation (∆H a ) and reaction (∆H r ) enthalpies, respectively.
In the upper panel, it can be seen that ∆H f has an almost linear behaviour with respect to n for both doped results, with a more pronounced slope in the case of hole doping. This 18 leads to more favourable formation enthalpies for hole doping at every n with respect to undoped results, which is particularly evident for n=5 (about 1.0 eV larger). A similar behaviour was found [57] in the case of di-hydrogenated (two hydrogen adatoms) graphene with hole doping leading to larger binding energies with respect to electron doping. In our calculation of the charge doped isolated circumcoronene it was found that, in both cases More interestingly, both electron and hole doping provide a reduction of ∆H a , which is particularly significant for n=4, 6, as seen in the intermediate panel.
As a result, for n > 3 the activation enthalpy ranges around 1.5 eV for both dopings, mimicking the behaviour reported for protons in Table 1 and Fig. 4 of Ref. [18], with electron doping in general leading to slightly lower values. In the case of the electron doping, we have checked for n=4 that the TS is characterized by an extra stabilization with respect to the reactant due to the larger affinity of the former. This is confirmed by a charge distribution analysis [58] ( see Table   S2 in the Supporting Material) which has shown that those carbon atoms of the TS supporting the chemisorbed hydrogen atoms and characterized by the sp 2 hybridization are indeed bearing a larger fraction of the negative charge with respect to the sp 3 hydridized neighboring atoms (as well as compared with their counterparts on the reactant species). As for the hole doping, both reactant and TS species are about 14-15 eV less stable than the corresponding neutral counterparts due to the high ionization energy of the latter. However, we have found (i.e. for n=4) a lower ionization energy in the case of the neutral TS of about 1.0 eV with respect to that of the neutral reactant and this can be due to its higher capacity of charge delocalization, as found by a larger partial positive charge globally bore by both the chemisorbed hydrogen and related supporting carbon atoms ( see Table S2).
On what concerns the effect on ∆H r , in the lower panel it can be noted that while electron doping is in general slightly more exothermic, the opposite can be observed for hole doping.
These results suggest that charge doping can help to reduce the barrier associated to the hydrogen flipping mechanism: in general, electron doping leads to slightly more favourable values for both activation and reaction enthalpies. In addition, we have checked that in the case of the electron doping, even if the activation enthalpy is also diminished for the diffusion and desorption processes, this reduction is less relevant and put the related barriers in the same range (around 1.4-1.7 eV) as for the flipping processes (see the results in Table   S3 of the Supporting Material obtained for the n=4 arrangement). This does not occur in the case of hole doping since this leads to a more significant decrease of the activation enthalpy (down to 0.9 eV) for the diffusion outcome which clearly becomes the most likely process (see Table S3). Therefore, we can conclude that when a given graphenic ring is

Isotopic effects
Since the reported activation and reaction enthalpies include ZPE contributions it is also possible to evaluate the corresponding thermodynamic properties for deuterium (D) chemisorbed atoms. To do that we have first considered the flipping of one deuterium atom for the n=4,5 cases (as in Fig.1 but with chemisorbed deuterium atoms ) and calculated the ∆H a and ∆H r enthalpy balances, which are compared with those for the non-deuterated counterparts, as shown in Table S4  is slightly favoured by just a few meV with respect to that of one deuterium atom. This small deviation is a consequence of the fact that ZPE corrections are important not only for the reactant but also for the transition and product states which provide comparable ZPE values, being their hydrogen/deuterium difference around 90 meV per chemisorbed atom, in agreement with previous determinations [50].
The obtained small differences demonstrate an almost negligible isotopic effect for the flipping process and they are certainly not sufficient to explain that deuterium permeation is not experimentally discerned in Ref. [20] within the considered detection limit. However, in addition to the flipping, we believe that other processes should be taken into account such as molecular hydrogen(deuterium) recombination. In fact, once at least two H(D) atoms found themselves on the other side of the graphene membrane they should react to form H 2 (D 2 ) in order to increase the container gas pressure which is responsible for the lifting of the membrane itself. Therefore, ZPE effects have been also considered for the recombination paths starting from the n=2 arrangements and through the para and ortho mechanisms (see Table III), which lead to H 2 (D 2 ) formation, and the related estimations are also reported in Table S4. In there it can be seen that H 2 formation is less activated than that of D 2 by about 65-100 meV and that the lighther isotope leads to a larger exothermic process of about 90-100 meV. We consider that the computed ∆H a and ∆H r estimations, which evidence non negligible differences for the recombination related to two isotopic species, clearly indicate that the permeation of the lighter species is globally favoured both kinetically and thermodynamically. We believe that present predictions on graphene isotopic selectivity are qualitatively compatible with the experimental findings of Ref. [20]. Moreover, additional quantum effects such as tunneling, here not considered for the sake of simplicity, are also expected to further facilitate [13] hydrogen transmission with respect to deuterium, therefore providing an additional support to the experimental evidences.

SUMMARY AND OUTLOOK
By means of DFT computations we have elucidated a cooperative mechanism leading to the flipping of a chemisorbed atomic hydrogen from one side to the other of a graphene layer.
We have demonstrated that the effectiveness of this mechanism, which involves a rearrangement of the orbital hybridization of the participant carbon atoms and the insertion of the flipping species between them, is highly dependent on the local density of the chemisorbed hydrogen atoms on the 2D carbon surface. In particular, for a low saturation of a given graphenic ring the transmission of atomic hydrogen is highly unlikely since alternative outcomes such as diffusion and desorption exhibit quite lower activation enthalpies. However, when a graphenic ring tends to saturation, the proposed flipping mechanism becomes competitive leading to a significant exothermic process. Moreover, we have shown that there exists at least one specific configuration of few neighboring chemisorbed hydrogen atoms for which the flipping of one hydrogen becomes the most likely process with a related limiting activation enthalpy which is in the range of the recent experimental estimation [20] of the barrier for hydrogen transmission through defect-free graphene. Present results provide solid support to a flipping mechanism responsible for the hydrogen transmission through pristine graphene and also suggest that electron doping could help to improve its performance.
We also consider that the proposed mechanism can be useful to provide a deeper understanding of atomic hydrogen intercalation which has been experimentally observed [24,25] in hydrogenated graphene/substrate interfaces. Additionally, ZPE effects have been also considered and they suggest that the trasmission of H 2 is favoured over D 2 both kinetically and thermodynamically. In order to fully assess the role of quantum effects on hydrogen trasmission in the near future we plan to investigate the impact of tunneling on the involved processes. We also believe that a detailed investigation of the effect of the membrane curvature on the computed magnitudes is also highly desirable and deserves to be addressed.  S1: Formation enthalpy (∆H f ) for the addition of n hydrogen atoms on a graphenic ring (see Eq. 1 in the main manuscript) as well as activation (∆H a ) and reaction (∆H r ) enthalpies for the flipping process as estimated by using the PBE, B3LYP and M062X functionals together with the cc-pVTZ basis sets. The n=4,5 arrangements (see Fig. 1    respectively. Σ H and Σ C are the sums of the charges of these four H and C atoms, respectively, and Σ H,C is the total sum. Labeling of the considered atoms corresponds to that in the image shown above, where the geometry of the transition state of the neutral prototype is depicted.

Reactants
Transition State q = −2 q = 0 q = +2 q = −2 q = 0 q = +2   : Isotopic dependence of the activation (∆H a ) and reaction (∆H r ) enthalpies for the flipping and recombination processes. For the flipping the n=4,5 initial arrangements of chemisorbed atoms (hydrogen or deuterium) have been considered while recombination refers to the n=2 case and occurring for both para and orho mechanisms (see also