# Chemisorption and physisorption of an H atom on a graphenic surface: Interaction potential and sticking probabilities

Mukul Dhiman, Sabine Morisset

► To cite this version:

Mukul Dhiman, Sabine Morisset. Chemisorption and physisorption of an H atom on a graphenic surface: Interaction potential and sticking probabilities. Chemical Physics Letters, 2025, 871, pp.142104. 10.1016/j.cplett.2025.142104. hal-05061300

HAL Id: hal-05061300

https://hal.science/hal-05061300v1

Submitted on 9 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Research paper

# Chemisorption and physisorption of an H atom on a graphenic surface: Interaction potential and sticking probabilities

Mukul Dhiman<sup>b,a</sup>, Sabine Morisset<sup>a,b</sup>,*

<figure><img src="image_3.png" /></figure>

## HIGHLIGHTS

* The energy dissipation (phonons) in the surface is included in the dynamics calculation.

* The H atom can stick in chemisorption and in physisorption on a graphenic surface.

* The sticking of H atom in chemisorption well can take place by tunnelling.

* The Van der Waals forces are included in the DFT calculations.

## GRAPHICAL ABSTRACT

<figure><img src="image_4.png" /></figure>

## ARTICLE INFO

**Keywords:**

Mixed classical-quantum dynamics

Wave packet propagation

Dissipation energy

Phonon

Sticking probabilities

## ABSTRACT

A complete study on the sticking of a hydrogen atom on the graphene surface has been performed, including energy dissipation in the surface via phonons. The sticking probabilities of H atom in physisorption and chemisorption wells are calculated for the first time in single dynamics calculation. A mixed classical-quantum dynamics method is used. These calculations are parametrized using a Van der Waals functional to obtain potential and phonon modes. In the astrophysical context, our results show that at a grain temperature of 10 K and small initial energy of the H atom, is predominately physisorbed in the v=1 vibrational state.

## 1. Introduction

The interaction of an H atom with the graphenic surface has been studied for many years. The single H atom sticking mechanism is a cornerstone for H₂ formation in the interstellar medium (ISM) [1], formed

following the Eley-Rideal mechanism and the Langmuir-Hinshelwood mechanism, as seen experimentally [2]. In the ISM, the interstellar dust grain can be in amorphous carbon, Polycyclic Aromatic Hydrocarbon (PAH, [3]) or cover of ices. H atoms sticking to dust grains can be

<footnote><sup>a</sup> Institut des Sciences Moléculaires d'Orsay, ISMO, UMR 8214, CNRS, Université Paris Saclay, F-91405 Orsay, France</footnote>

<footnote><sup>b</sup> MSME, Université Gustave Eiffel, UPEC, CNRS, UMR 8208, F-77454 Marne-La-Vallée, France</footnote>

<footnote>* Corresponding author.
E-mail address: sabine.morisset@cnrs.fr (S. Morisset).</footnote>
a crucial means of cooling interstellar gas [4]. This mechanism is also essential in hydrogen storage devices, where interaction between functionalized graphene with an H atom or H₂ molecule has been studied [5–7].

In these domains of interest, our study is focused on the sticking of one H atom on a single layer of graphene. Multilayer graphene seems to reduce the binding energy [8].

In this work we look at the H atom approaching one of the carbon atom collinearly. The potential of such system shows an existence of a long range physisorption well followed by a barrier and a short range chemisorption with surface reconstruction — the collinear C atom protrudes through the surface to create a chemical bond with the H atom [9–13].

Energy dissipation through phonon modes is essential to study the adsorption dynamics [14–17]. The sticking dynamics of the H atom interacting with a graphenic surface has been studied for many years. Phonon mediated sticking dynamics in chemisorption has been studied using close coupling wave packet method [18], mixed classical-quantum method in the one-phonon exchange approximation [19,20] and classical, quasi-classical and quantum mechanical methods [21]. Cazaux et al. [20] have shown that at small energy the sticking probability is low and is dependent on the surface temperature. Bonfanti et al. [21], at T = 0 K, show that a large part of the energy is transferred to the phonon bath, which implies that the multi-phonon process is operative, especially for the low frequencies. In physisorption, Lepetit et al. [22,23] have studied quantum mechanically the sticking of H atom on supported and suspended graphene surface. They conclude that at 10 K, the sticking is efficient at energy smaller than 10 meV and that the sticking probability decreases with increasing energy. The sticking of H atom on grains is studied in either chemisorption or physisorption separately, but not together. The objective of this work is to study, in the same dynamics calculation, the sticking of an H atom on graphenic surface, including the chemisorption and the physisorption wells.

The first part describes the H atom interaction with the surface using the Density Functional Theory (DFT) including the non-local van der Waals (vdW) functional for long-range interactions. The choice of the functional is discussed, as well as phonon calculations at the gamma point in the Irreducible Brillouin Zone (IBZ). The second part is devoted to the dynamics calculations using the mixed classical-quantum method [19] - classical phonons and a wave packet propagation method for the motion of H atom.

## 2. Potential and phonon calculations

### 2.1. Potential calculation

The H atom interaction with the graphenic surface is calculated using the DFT method. These DFT calculations are performed using the Plane Wave Self-Consistent Field (PWSCF) method using the Quantum Espresso (QE) package [24–26]. In the PWSCF code, the non-local term is added using the vdW functional for the description of the exchange-correlation terms. Ultra-soft pseudopotentials are used [27]. The optimized kinetic energy cut-off is 100 Ry for the wave function and 800 Ry for the charge density for vdW SCF calculations. The spin polarization scheme is used. Calculations are carried out for a 2 × 2 graphene supercell containing 8 carbon atoms and one H atom, using 8 × 8 × 1 k-points. The energy convergence is less than 10meV. The total energy $E_{int}$ of the super-cell is calculated as follows:

$$E_{int} = E_{H-graph} - E_{graph} - \frac{1}{2} E_{H_2} \quad (1)$$

where $E_{H-graph}$ is the energy of interaction of H atom with the surface, $E_{graph}$ is the graphenic surface alone and $E_{H_2}$ is the energy of H₂ molecule.

Previous calculation show that the H atom can be physisorbed at 3.0 Å from the graphene surface [9,10] and at small distance between

<figure><img src="image_5.png" /><figcaption>Fig. 1. Potential Energy Curve for co-linear H atom on graphene using varied DFT functional on 2 × 2 graphene supercell. The insert shows the physisorption plot.</figcaption></figure>

1.45–1.50 Å it is chemisorbed [12,13,28,29]. Using MP2, Bonfanti et al. [10] obtain a physisorption well of 34.5 meV, which is in good agreement with the experimental values of 39.2 ± 0.5 meV of Ghio et al. [9]. The chemisorbed well is about 0.76 eV, a chemical bond is created with the C atom and a surface reconstruction is observed, where the C atom goes out of the surface plane: hybridization sp2 to sp3 [12,13,28,29].

In this work, the interaction potential is calculated using the PBE functional, and we have tested several forms of the non-local term added to the PBE functional: vdW-DF [30,31], vdW-DF-CX [32], and vdW-RVV10 [33] non-local functionals. Non-local vdW density functionals are defined in terms of the electron density alone, without any orbital input. These non-local functionals are very general, as they require neither splitting the system into interacting fragments, nor do they require any kind of atomic partitioning. They treat inter- and intra-molecular dispersion interactions on equal footing. Berland et al. [32] (vdW-DF-CX) improve upon the vdW-DF given by Dion et al. [30,31] by enhancing the exchange accuracy, while vdW-RVV10 is the revised VV10 [34] pseudopotential given by Sabatini et al. [33], which further improves computational efficiency and accuracy in describing vdW interactions.

Fig. 1 and Table 1 present a comparison of different interaction potentials obtained using PBE alone and with vdW forces. Using the PBE functional, the physisorption well is not obtained, and the chemisorption well (−0.784 eV) and barrier (0.18 eV) are in good agreement with previous studies [8,12]. Using the PW91 functional [12], a small value for the physisorption well is obtained (−31 meV), whereas with the PBE functional, the physisorption well is not obtained. The PBE and PW91 functionals are known to exhibit a large discrepancy in results when calculating properties where surface effects are present [35]. As the PW91 functional does not include the non-local term for long-range interactions, the value obtained with this functional may vary at large distances. Using different vdW functionals, the physisorption well at large distances is shown in Fig. 1. Using non-local DFT functionals, see Table 1, though we recover the long-range properties, we underestimate the chemisorption well and the adsorption barrier height. Following that, the best results are obtained with the vdW-RVV10 functional ($E_{barrier}$ = 0.169 eV, $E_{chem}$ = −0.765 eV, and $E_{phys}$ = −0.0294 eV), which were also reported by Bonfanti et al. [8]. In the two calculations done independently, the physisorption well is correctly represented with the vdW-RVV10 functional compared to previous works [12,36].

### 2.2. Phonon calculations

The vibrational modes (phonons) of the H atom interacting with the graphene surface are calculated using the phonon module in the
<caption>**Table 1**</caption>

Barrier Height (E<sub>barrier</sub>), Chemisorption Energy (E<sub>chem</sub>) and Physisorption Energy (E<sub>phys</sub>) in eV as calculated using DFT in Quantum Espresso 7.2 for a single Hydrogen atom on Graphene surface.

<table><thead><tr><th rowspan="2">Functional</th><th colspan="2">E<sub>barrier</sub> (eV)</th><th colspan="2">E<sub>chem</sub> (eV)</th><th colspan="2">E<sub>phys</sub> (eV)</th></tr><tr><th>This work</th><th>Lit</th><th>This work</th><th>Lit</th><th>This work</th><th>Lit.</th></tr></thead><tbody><tr><td>PW91 [37]</td><td>-</td><td>0.185 [12]</td><td>-</td><td>-0.608 [12]</td><td>-</td><td>-0.031[12]</td></tr><tr><td>PBE [38]</td><td><b>0.188</b></td><td>0.3 [8]</td><td><b>-0.784</b></td><td>-0.8 [8]</td><td>No</td><td>No [8]</td></tr><tr><td>vdW-RVV10 [33]</td><td><b>0.169</b></td><td>0.179 [8]</td><td><b>-0.765</b></td><td>-0.851 [8]</td><td><b>-0.0294</b></td><td>-0.029 [8]</td></tr><tr><td>vdW-DF-CX [32]</td><td><b>0.031</b></td><td>-</td><td><b>-0.905</b></td><td>-</td><td><b>-0.0677</b></td><td>-</td></tr><tr><td>vdW-DF [30,31]</td><td><b>0.107</b></td><td>-</td><td><b>-0.812</b></td><td>-1.10 [36]</td><td><b>-0.0645</b></td><td>-0.020 [36]</td></tr><tr><td>MP2</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-0.0345 [10]</td></tr></tbody></table>

Quantum Espresso (QE) package [24–26]. The objective is to calculate the phonon modes using the vdW functionals and compare them with previous results [18] obtained using PBE functional. Firstly, the phonon are modes calculated at the gamma point in the IBZ of a single graphene sheet which are decomposed into 3 acoustic modes identified as: in-plane longitudinal modes (LA), in-plane transverse modes (TA) and out-of-plane transverse modes (ZA). Similarly, there are 3 optical phonon modes of graphene are present around 1550 cm⁻¹ (LO), 1550 cm⁻¹ (TO) and 850 cm⁻¹ (ZO). Following that, the perturbed phonon modes are obtained where the H atom is interacting with the surface. The three modes associated with the chemisorbed H atom on the graphene are: the first is around 2600 cm⁻¹, which is the associated with out-of-plane vibration of H atom on the graphene and the other two degenerate modes around 1200 cm⁻¹, which are the vibration of H in the plane parallel to graphene. Using the PBE functional, all phonon modes of the H atom chemisorbed on the graphene are reported in Fig. 2. These results are in reasonable agreement with experiments of Zecho et al. [39], Allouche et al. [29] and previous calculations of Morisset et al. [18].

Table 2 gives the phonon modes obtained using different vdW functionals in comparison with phonons obtained using PBE functional. Several phonon modes are not converged for most of the vdW functionals. The best result of all is obtained using the vdW-RVV10 functional, where the lowest 3 phonon modes can be neglected for periodic systems and other frequencies are within 50 cm⁻¹ of PBE phonon modes. The non-converged frequencies are removed by applying the crystal Acoustic Sum Rule (ASR). The phonon frequencies are obtained using the vdW-RVV10 functional at gamma point in the IBZ. Fig. 2 shows the overlay of the potential and associated phonon frequencies after application of ASR. These frequencies are used in dynamics calculations.

The remaining negative frequencies observed in Fig. 2 could point to some physical phenomenon. In the following case, at the chemisorption barrier, there are negative frequencies of -430 cm⁻¹ and -281 cm⁻¹ which point to the presence of the transitions states. At large distances, the vibrational modes are those of the graphene alone. The chemisorption well has been previously noted and is located at 1.5 Å from the graphene sheet with degenerated in plane frequencies of 1138 cm⁻¹ and 2616 cm⁻¹ out of plane associated to C-H vibration. The chemisorption barrier is located at around 2 Å with degenerated in frequencies of 260 cm⁻¹ in plane and one -430 cm⁻¹ out of plane associated to the C-H bonding transition state. In this work, we also observe a physisorption well located at 3 Å and the frequencies of the H atom are 254 cm⁻¹ out of the plane and degenerate 224 cm⁻¹ in the plane. In all three modes, the carbon atom is observed to be affected by the presence of the H atom. For small z values, the high frequency of about 3000 cm⁻¹ is observed, which can be associated with high vibrational excitation of C-H. We also suspect the presence of a low-lying transition state at a large distance of about 3.5 Å where we have a -281 cm⁻¹ mode.

<caption>**Table 2**</caption>

Phonons frequencies calculated on 2 × 2 working cell for the H atom chemisorbed on the single layer graphene using the phonon code of Quantum Espresso package [24–26] obtained as a function of the functional.

<table><thead><tr><th>number</th><th>vdW-DF</th><th>vdW-DF-CX</th><th>vdW-RVV10</th><th>PBE</th></tr></thead><tbody><tr><td>1</td><td>-864.900109</td><td>-774.443572</td><td>-155.843572</td><td>49.213912</td></tr><tr><td>2</td><td>-864.900109</td><td>-774.443537</td><td>-102.064976</td><td>49.213912</td></tr><tr><td>3</td><td>-745.741165</td><td>-616.43166</td><td>-102.064828</td><td>103.898705</td></tr><tr><td>4</td><td>-676.973448</td><td>-275.761449</td><td>343.953335</td><td>395.159297</td></tr><tr><td>5</td><td>-644.335961</td><td>-275.760561</td><td>343.953623</td><td>395.170684</td></tr><tr><td>6</td><td>-644.335961</td><td>168.061103</td><td>449.67022</td><td>461.286595</td></tr><tr><td>7</td><td>-439.485144</td><td>418.206046</td><td>564.917395</td><td>590.724811</td></tr><tr><td>8</td><td>-439.485144</td><td>418.206102</td><td>564.917478</td><td>590.728089</td></tr><tr><td>9</td><td>416.831269</td><td>470.125596</td><td>591.6489</td><td>601.297983</td></tr><tr><td>10</td><td>464.284689</td><td>629.425177</td><td>611.885963</td><td>621.656535</td></tr><tr><td>11</td><td>464.284689</td><td>768.41307</td><td>612.120676</td><td>621.663507</td></tr><tr><td>12</td><td>541.173071</td><td>768.413107</td><td>612.120716</td><td>626.196165</td></tr><tr><td>13</td><td>674.070333</td><td>799.12411</td><td>839.760371</td><td>855.185764</td></tr><tr><td>14</td><td>788.313219</td><td>966.88121</td><td>1125.08559</td><td>1127.003002</td></tr><tr><td>15</td><td>931.142905</td><td>966.88128</td><td>1125.085684</td><td>1127.011600</td></tr><tr><td>16</td><td>931.142905</td><td>989.872528</td><td>1239.859871</td><td>1248.649052</td></tr><tr><td>17</td><td>1147.185712</td><td>1154.133163</td><td>1239.859997</td><td>1248.672517</td></tr><tr><td>18</td><td>1196.487048</td><td>1222.888953</td><td>1263.246687</td><td>1259.547605</td></tr><tr><td>19</td><td>1196.487048</td><td>1222.889074</td><td>1308.575435</td><td>1321.394910</td></tr><tr><td>20</td><td>1211.813237</td><td>1226.919149</td><td>1321.336328</td><td>1336.777999</td></tr><tr><td>21</td><td>1211.813237</td><td>1226.919561</td><td>1337.235985</td><td>1336.782407</td></tr><tr><td>22</td><td>1261.64491</td><td>1265.714862</td><td>1337.236255</td><td>1339.981238</td></tr><tr><td>23</td><td>1328.227509</td><td>1333.27116</td><td>1376.559323</td><td>1388.931094</td></tr><tr><td>24</td><td>1328.227509</td><td>1333.271409</td><td>1376.559731</td><td>1388.940513</td></tr><tr><td>25</td><td>1394.508292</td><td>1403.712742</td><td>1524.152361</td><td>1534.732433</td></tr><tr><td>26</td><td>1394.508292</td><td>1403.713072</td><td>1524.152584</td><td>1534.750737</td></tr><tr><td>27</td><td>2566.231309</td><td>2568.204887</td><td>2611.877555</td><td>2628.310468</td></tr></tbody></table>

### 3. Dynamics calculations

#### 3.1. Brief description of the dynamics method

The detailed theory is presented in the previous papers of Morisset et al. [18–20]. The Hamiltonian used to describe the H atom in interaction with the dynamical surface in the one-phonon approximation exchange is:

$$H = T + V_{\text{tot}} + H_{\text{bath}} \qquad (2)$$

where T is the kinetic operator, V<sub>tot</sub> is the gas-surface potential operator [18–20] and H<sub>bath</sub> is the Hamiltonian for the lattice vibration without coupling to an atom:

$$H_{\text{bath}} = \sum_{i=1}^{N} \frac{1}{2} \omega_i(z) (Q_i^2 + P_i^2) \qquad (3)$$

where N is the number of phonon modes Q<sub>i</sub> and P<sub>i</sub> are the position and momenta of the vibrational modes of frequency ω<sub>i</sub>. V<sub>tot</sub> is the sum of the static potential describing the collinear approach of the H atom to the graphenic surface V<sub>0</sub>(z) and the interaction potential between the H atom motion coupled to the phonon bath:

$$V_{\text{tot}} = V_0(z) + \sum_{i=1}^{N} g_i(z) Q_i \qquad (4)$$

where g<sub>i</sub>(z) is the coupling term between the H motion and the phonon bath [18]. The coupling term, g<sub>i</sub>(z) on the z grid, is calculated by interpolating the dynamical matrix obtained by the phonon code of QE package [24–26].
<table><caption>**Table 3**
Sticking probabilities as a function of surface temperature T and of the initial energy E₀ of the wave packet, N<sub>traj</sub> = 400.</caption><thead><tr><th rowspan="2">E<sub>0</sub> (meV)</th><th colspan="3">50</th><th colspan="3">177</th></tr><tr><th>10</th><th>100</th><th>300</th><th>10</th><th>100</th><th>300</th></tr></thead><tbody><tr><th>P<sub>S</sub></th><td>1.8500×10<sup>-2</sup></td><td>0.07770</td><td>0.1923</td><td>0.333</td><td>0.461</td><td>0.584</td></tr><tr><th>P<sub>Chem</sub></th><td>0.0026×10<sup>-2</sup></td><td>0.00521</td><td>0.0341</td><td>0.216</td><td>0.338</td><td>0.448</td></tr><tr><th>P<sub>Phys</sub></th><td>1.8470×10<sup>-2</sup></td><td>0.07250</td><td>0.1581</td><td>0.117</td><td>0.128</td><td>0.135</td></tr></tbody></table>

<table><caption>**Table 4**
Percentage of the sticking probability on each vibrational states in chemisorption and physisorption.</caption><thead><tr><th rowspan="2"></th><th colspan="4">Chemisorption</th><th colspan="3">Physisorption</th></tr><tr><th>v = 0</th><th>v = 1</th><th>v = 2</th><th>v = 3</th><th>v = 0</th><th>v = 1</th><th>v = 2</th></tr></thead><tbody><tr><th>54 meV 10 K</th><td>&lt;0.01%</td><td>&lt;0.01%</td><td>&lt;0.01%</td><td>0.13%</td><td>5.1%</td><td>92.9%</td><td>1.5%</td></tr><tr><th>54 meV 100 K</th><td>2.6%</td><td>0.8%</td><td>1.6%</td><td>1.6%</td><td>14.5%</td><td>73.6%</td><td>5.0%</td></tr><tr><th>54 meV 300 K</th><td>5.3%</td><td>1.6%</td><td>5.0%</td><td>5.7%</td><td>22.3%</td><td>55.2%</td><td>4.8%</td></tr><tr><th>177 meV 10 K</th><td>22.8%</td><td>12.0%</td><td>14.0%</td><td>16.2%</td><td>13.1%</td><td>7.7%</td><td>14.3%</td></tr><tr><th>177 meV 100 K</th><td>22.5%</td><td>23.1%</td><td>11.7%</td><td>15.8%</td><td>10.6%</td><td>5.9%</td><td>10.0%</td></tr><tr><th>177 meV 300 K</th><td>20.9%</td><td>24.7%</td><td>15.0%</td><td>16.2%</td><td>11.7%</td><td>4.7%</td><td>6.8%</td></tr></tbody></table>

The sticking of the H atom on the surface is treated using a mixed quantum-classical method. The motion of the H atom along the z-coordinate is treated quantum mechanically using the wave-packet propagation method and each vibrational mode is treated classically by solving the equations of motion derived from the Hamilton equation:

$$ \dot{Q}_i = \frac{\partial H_{\text{eff}}}{\partial P_i} \qquad (5) $$

$$ \dot{P}_i = -\frac{\partial H_{\text{eff}}}{\partial Q_i} \qquad (6) $$

where $H_{\text{eff}}$ is the effective Hamiltonian, given by the expectation value of the Hamiltonian $H$:

$$ H_{\text{eff}} = \langle \psi | H | \psi \rangle \qquad (7) $$

where $\psi$ is the wave function of the H atom along the z coordinate. At $t=0$, $\psi$ is a Gaussian wave packet with initial energy $E_0$, detailed in previous work [18]. Lanczos method is used for the propagation of the wave packet [40]. The kinetic energy operator is evaluated using the Fourier transform method [41]. All vibrational modes obey the equipartition theorem, that is, in each mode there is $\frac{k_B T}{2}$ of energy, where $T$ is the surface temperature. The initial temperature distribution is randomly selected, as suggested by Billing [42]. N<sub>traj</sub> gives the number of initial trajectories to reach convergence [19]. At the end of the propagation, the wave function is analysed using the reactive flux for the calculation of the total sticking probability ($P_S$) as described in previous papers [19,20].

## 3.2. Sticking probabilities

The sticking probabilities are calculated for the mean energies of the initial wave packet $E_0$ of 54 meV and 117 meV for temperatures of 10 K, 100 K and 300 K. In Tables 3 and 4, the sticking probabilities and their distribution over vibrational modes are reported. These results are converged using 400 trajectories and a dynamics time of 1.5 ps. Probabilities in chemisorption ($P_{Chem}$) and physisorption ($P_{Phys}$) potential wells are deduced from the projection of the wave function on the vibrational state of the chemisorption and physisorption wells. To obtain the probability in each state, this percentage should be multiplied by $P_S$.

For energy, $E_0 = 54$ meV, smaller than the value of the chemisorption barrier, $P_S$ increases with the temperature. As shown in Table 3, $P_{Chem}$ increases with the temperature and is smaller than $P_{Phys}$. In Table 4, for $T = 10$ K, the percentage of $P_S$ to find the H atom in vibrational state $v = 1$ in the physisorption well is about 93%. With increasing

temperature, the shared probabilities with other vibrational states ($v = 0$ and 2) are increased. However, the probability of observing the H atom in chemisorption remains low. Hence, at small energy and low temperature (10 K), the H atom is most likely to be physisorbed in the $v = 1$ vibrational state. Therefore, for collision energies smaller than the value of the chemisorption barrier, the H atom is preferably physisorbed on the graphenic surface, independent of the temperature. In Fig. 3 (bottom), at $T = 300$ K and $E_0 = 54$ meV, the sticking probabilities ($P_S$, $P_{Chem}$ and $P_{Phys}$) as function of time and temperature are reported where the state $v = 1$ of the physisorption well is mostly populated. At 54 meV, the tunnelling between physisorption and chemisorption wells is not efficient.

For energy, $E_0 = 177$ meV, higher than the energy of the chemisorption barrier, $P_{Chem}$ and $P_S$ increase with the temperature while, $P_{Phys}$ almost remains constant (Tables 3 and 4). The sticking of the H atom is dominated by the sticking in chemisorption at this energy. But the sticking in physisorption is not negligible. Moreover, the results presented in this work show that the sticking probabilities increase with the temperature, whatever the initial energy of the H atom. The sticking probability in the physisorption is observed to be almost equally distributed over the three vibrational modes and is independent of the temperature. In the chemisorption well, we observe almost the same behaviour. At $T = 10$ K, the percentage of the sticking probability in state $v = 0$ is highest, and with increasing temperature, the sticking probability in the state $v = 1$ increased while in other states it did not change much in chemisorption, while the sticking probability in physisorption is reduced. Hence, as expected, with increasing temperature, at higher energy, the sticking probability in physisorption is reducing while in the chemisorption it is increasing. In Fig. 3 (top), the sticking probabilities as a function of time and temperature are reported for 100 K and 177 meV. We see that $P_S$ does not change with time, but the probability of chemisorbed H atom increases while physisorbed H atom reduces. This could be due to tunnelling between the two potential wells.

In previous work of Cazaux et al. [20], the sticking probability in chemisorption at 10 K noted to be about 0.8 at energy about 170 meV, and about $3 \times 10^{-4}$ at energy about 60 meV. Similar to previous work [20], this work uses the identical approach for the dynamics calculations and the model used to calculate couplings between the phonon modes and motion of the H atom. Hence, the sticking probabilities can be directly compared. In Table 3, we see that the sticking probability at 54 meV is $1.85 \times 10^{-2}$, which is higher than previous results [20]. This increase of $P_S$ is due to the physisorbed H atom,

<figure><img src="image_6.png" /><figcaption>**Fig. 2.** In red (left y-axis) phonon frequencies at gamma point in the Irreducible Brillouin Zone as a function of the distance of H atom to the surface using vdW-RVV10 functional after application of Acoustic Sum Rule. In blue (right y-axis), the interaction potential with vdW-RVV10 functional as a function of the distance of the H atom to the surface. The barrier is at 2 Å, the chemisorption well is at 1.5 Å and the physisorption well is at 3 Å.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Fig. 3.** Top: Sticking probabilities at T = 100 K and initial energy of 177 meV as function of the time propagation. Bottom: Sticking probabilities at T = 300 K and initial energy of 54 meV as function of the time propagation. The blue line with hollow squares is the total sticking probability P<sub>S</sub>. The red line with hollow circles is the sticking probability in chemisorption well P<sub>Chem</sub> and Black line with the cross is the sticking probability in physisorption well P<sub>Phys</sub>. The left green square, diamond, triangle and circle are the sticking probabilities in state v = 0, v = 1, v = 2, v = 3, respectively in the chemisorption well and right are the sticking probabilities in state v = 0, v = 1, v = 2 respectively in the physisorption well.</figcaption></figure>

which is not observed in previous work. However, for energy greater
than energy of chemisorption barrier, 177 meV, P<sub>S</sub> is 0.333 at 10 K
and 0.461 at 100 K. These sticking probabilities are smaller than the
previous results [20]. These discrepancies can be due to following
reasons: (i) The H atom is physisorbed on the surface, and it stays
in this state, implying inefficient or negligible tunnelling between the
chemisorption and the physisorption potential wells (Table 4); (ii) The
other difference could be due to the missing vibrational modes in the
IBZ. In this work, only gamma point vibrational modes are included
in the IBZ, while in previous work [20] 61 vibrational modes in the
IBZ were included. The excess of vibrational modes in IBZ could aid in
increasing the possibility of energy exchange between the incident H
atom and phonon modes. Perhaps, the sticking probabilities reported
in this work are underestimated. However, using the limited phonon
calculations on the gamma point in IBZ, this work still captures the H
sticking dynamics in the chemisorption and physisorption well on the
graphene layer.

Lepetit et al. [23] have obtained, using dynamics calculation, a
sticking probability of a few percent for collision energy between 0
and 25 meV, increasing with the temperature. As discussed previously,
at small incident energy, the H atom preferentially sticks at a large
distance on the graphenic surface. The sticking probability at these
small energies has the same behaviour as the one calculated with
only the physisorption well, shown in the work of Lepetit et al. [23].
Moreover, they [23] have also reported the individual sticking prob-
abilities at T = 10 K and T = 300 K in the energy domain of 0 to
25 meV. They show peaks of resonances which illustrate the sticking
and trapping in different vibrational states of the stuck H atom, per-
pendicular and parallel to the surface. They reported the states, v = 0
and v = 1 of the stuck H atoms perpendicular to the surface are the
most important for sticking probabilities, and the sticking probabilities
are independent of the temperature. Although it is challenging to
make a direct comparison between Lepetit et al. [23] and our results
due to differing methodologies, particularly the inclusion of surface

corrugation, our results also demonstrate that the H atom is stuck in the
lowest vibrational states of the chemisorption and physisorption wells.

Experimentally, Zecho et al. [39] used an atomic beam with energy
centred at about 0.2 eV, and they obtained a sticking coefficient of
about 0.5 for a saturation coverage of around 0.3 with a surface
temperature of 150 K. In this experiment, the sticking probabilities are
dependent on the saturation coverage (0.4 ± 0.2) and vary between
0.25 and 0.5. Although the results presented in this work are difficult
to compare directly with the work of Zecho et al. [39] given that the
chemisorption barrier height is underestimated in this work, at 0.17 eV.
The nearest calculations with the experiment are those with the main
energy of 0.177 eV at 100 K, where the sticking probability is about
0.48. This result is in the same domain of results as Zecho et al. [39].
Our results seem to be good considering, only one-phonon processes
have been taken into account in this model and the cell size for the DFT
calculation is small, the coverage is of $\frac{1}{8}$ - there are 8 carbon atoms
and 1 hydrogen atom. Another experimental work [43] observed a
physisorption channel on the graphene surface from which the H atom
can be chemisorbed. In these dynamics, tunnelling from physisorbed H
atom to chemisorbed H atom is observed, see Fig. 3. The effect observed
is weak in this calculation as only a single phonon exchange is taken
into account.

4. Conclusion

We present for the first time in single dynamics calculation the
sticking probability of the H atom in physisorption and chemisorp-
tion wells. These mixed quantum classical dynamics calculations are
parametrized using DFT-vdW-RVV10 functional to obtain potential and
phonon modes.

Several vdW functionals have been tested and with the vdW-RVV10
functional - the barrier, the chemisorption and the physisorption en-
ergies are in good agreement with previous works [8,11,12,28,29,
39].
For small energy (55 meV), the H atom is physisorbed for all temperatures. The probability distribution of the physisorbed H atom at low temperature is on the $v = 1$ vibrational state. With increasing temperature, the distribution is shared among other vibrational states, but it is majorly located on the two lowest vibrational states. We also note the small, but non-negligible increase of the chemisorbed H atom with increasing temperature at low energy. At energy higher than the value of the chemisorption barrier (177 meV), the H atom is chemisorbed for all temperatures with a non-negligible probability for it to be physisorbed. The probability distribution at low temperature is highest for the lowest vibrational state $v = 0$. With increasing temperature, the probability is majorly distributed among the two lowest vibrational states. The distribution in the higher vibrational states of the chemisorption well increases marginally with temperature, while in all the vibrational states of the physisorption well it reduces marginally. Hence, we can say that with increasing temperature, the physisorbed H atom's sticking probability reduces while the chemisorbed H atom's probability increases, independent of the initial energy. We also report the visible tunnelling effect observed at 177 meV. We expect these effects to be better reproduced by incorporating multi-phonon exchange in the dynamics.

In the astrophysical context, our results show that at grain temperature of 10 K and small initial energy of the H atom, is predominately physisorbed in the $v = 1$ vibrational state. With increasing initial energy, the sticking probability moves to chemisorption well. In this case, the sticking wave function is a combination of the two lowest vibrational states. These results are very important for the H₂ formation as seen in the Kinetic Monte Carlo simulation used in the previous paper of Cazaux et al. [20]. The study of sticking dynamics at room temperature (300 K) is relevant for hydrogen storage devices, where we show that a single H atom can be either physisorbed or chemisorbed on the graphene surface. An experimental study [43] has shown that, at high temperatures and low initial energy, the H atom will have enough thermal energy to be chemisorbed. This work confirms their observations using mixed quantum-classical dynamics. It is also known that the potential around changes after chemisorbed H atom, which facilitates the sticking of another H atom which might be physisorbed at another site(s), leading to the formation of H₂ molecule [20]. An interesting study following this would be to simulate the desorption dynamics of physisorbed and chemisorbed H atom. Also, using this model, it would be interesting to investigate the sticking probabilities in a multi-layer graphene structure.

In our approximation, we have only included the phonon modes at the gamma point in the IBZ and only single phonon exchange process for dynamics. In this case, the energy exchange between the surface and the H atom is lower than in previous calculations [20], and thus sticking probabilities are underestimated, especially in the chemisorption case where multi-phonon process are operative, especially at low frequencies [21]. Despite this approximation, our calculations demonstrate that the sticking in chemisorption and physisorption is important as a function of the initial energy of the H atom and the surface temperature.

### CRediT authorship contribution statement

**Mukul Dhiman:** Writing – original draft, Methodology, Formal analysis. **Sabine Morisset:** Writing – original draft, Supervision, Methodology, Formal analysis.

### Declaration of competing interest

The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: Morisset reports financial support was provided by French National Research Agency. If there are other authors, they declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgement

The authors acknowledge the support of the French Agence Nationale de la Recherche (ANR), under grant ANR-19-CE30-0039 (project HyTRAJ).

### Data availability

Data will be made available on request.

### References

[1] R. Gould, E. Salpeter, Astrophys. J. 138 (1963) 393.

[2] K.W. Kolasinski, Surf. Sci.:Foundations of Catalysis and Nanoscience, Wiley, Hoboken, N.J., 2008.

[3] J. Weingartner, B.T. Draine, Astrophys. J. 548 (2001) 296.

[4] M. Spaans, J. Silk, Astrophys. J. 538 (2000) 115.

[5] K. Iyakutti, R.P. Reji, K. Ajaijawahar, R.R. I. Lakshmi, V.J. Surya, A. Karthigeyan, Y. Kawazoe, Int. J. Quantum Chem. 124 (2024) e27467.

[6] V. Jain, B. Kandasubramanian, J. Mater. Sci. 55 (2020) 1865.

[7] V. Tozzini, V. Pellegrini, Phys. Chem. Chem. Phys. 18 (2013) 80.

[8] M. Bonfanti, R. Martinazzo, Phys. Rev. B 97 (2018) 117401.

[9] E. Ghio, L. Mattera, C. Salvo, F. Tommasini, U. Valbusa, Vibrational spectrum of H and D on the (0001) graphite surface from scattering experiments, J. Chem. Phys. 73 (1) (1980) 556-561.

[10] M. Bonfanti, R. Martinazzo, G.F. Tantardini, A. Ponti, J. Phys. Chem. C 111 (16) (2007) 5825-5829.

[11] N. Rougeau, D. Teillet-Billy, V. Sidis, Phys. Chem. Chem. Phys. 13 (2011) 17579-17587.

[12] Y. Ferro, F. Marinelli, A. Allouche, J. Chem. Phys. 116 (18) (2002) 8124-8131.

[13] X. Sha, B. Jackson, Surf. Sci. 496 (2002) 318.

[14] B. Smits, M. Sah, K. Naskar, S. Adhikari, J. Meyer, M.F. Somers, J. Chem. Phys. (161) (2024) 154702.

[15] M. Sah, K. Naskar, S. Adhikari, B. Smits, J. Meyer, M.F. Somers, J. Chem. Phys. (161) (2024) 014306.

[16] B. Jackson, J. Chem. Phys. (97) (1992) 6792.

[17] B. Jackson, J. Chem. Phys. (88) (1987) 1383.

[18] S. Morisset, A. Allouche, J. Chem. Phys. 129 (2008) 024509.

[19] S. Morisset, Y. Ferro, A. Allouche, J. Chem. Phys. 133 (2010) 044508.

[20] S. Cazaux, S. Morisset, M. Spaans, A. Allouche, Astron. Astrophys. A27 (2011) 10.

[21] M. Bonfanti, B. Jackson, K. Hughes, I. Burghard, R. Martinazzo, J. Chem. Phys. 143 (2015) 124704.

[22] B. Lepetit, B. Jackson, Phys. Rev. Lett. 107 (2011) 236102.

[23] B. Lepetit, D. Lemoine, Z. Medina, B. Jackson, J. Chem. Phys. 134 (2011) 114705.

[24] 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. Fabris, G. Fratesi, S. de Gironcoli, 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, R.M. Wentzcovitch, J. Phys.: Condens. Matter. 21 (2009) 395502.

[25] P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M.B. Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo, A.D. Corso, S. de Gironcoli, P. Delugas, R.A. DiStasio Jr., A. Ferretti, A. Floris, G. Fratesi, G. Fugallo, R. Gebauer, U. Gerstmann, F. Giustino, T. Gorni, J. Jia, M. Kawamura, H.-Y. Ko, A. Kokalj, E. Küçükbentli, M. Lazzeri, M. Marsili, N. Marzari, F. Mauri, N.L. Nguyen, H.-V. Nguyen, A.O. de-la Roza, L. Paulatto, S. Poncé, D. Rocca, R. Sabatini, B. Santra, M. Schlipf, A.P. Seitsonen, A. Smogunov, I. Timrov, T. Thonhauser, P. Umari, N. Vast, X. Wu, S. Baroni, J. Phys.: Condens. Matter. 29 (2017) 465901.

[26] https://www.quantum-espresso.org/.

[27] A.M. Rappe, K.M. Rabe, E. Kaxiras, J.D. Joannopoulos, Phys. Rev. B 41 (1990) 1227-1230.

[28] L. Jeloica, V. Sidis, Chem. Phys. Lett. 300 (1999) 157.

[29] A. Allouche, Y. Ferro, T. Angot, C. Thomas, J. Layet, J. Chem. Phys. 123 (2005) 124701.

[30] M. Dion, H. Rydberg, E. Schröder, D.C. Langreth, B.I. Lundqvist, Phys. Rev. Lett. 92 (2004) 246401.

[31] T. Thonhauser, S. Zuluaga, C.A. Arter, K. Berland, E. Schröder, P. Hyldgaard, Phys. Rev. Lett. 115 (2015) 136402.

[32] K. Berland, P. Hyldgaard, Phys. Rev. B 89 (2014) 035412.

[33] R. Sabatini, T. Gorni, S. de Gironcoli, Phys. Rev. B 87 (2013) 041108.

[34] O.A. Vydrov, T. Van Voorhis, J. Chem. Phys. 133 (24) (2010) 244103.

[35] A. Mattsson, R. Armiento, P. Schultz, T. Mattsson, Phys. Rev. B 73 (2006) 195123.
[36] M. Moaied, J.A. Moreno, M.J. Caturla, F. Ynduráin, J.J. Palacios, Phys. Rev. B 91 (2015) 155419.

[37] J. Perdew, J. Chevary, S. Vosko, K. Jackson, M. Pederson, D. Singh, C. Fiolhais, Phys. Rev. B 23 (1992) 6671.

[38] J.P. Perdew, M. Ernzerhof, K. Burke, J. Chem. Phys. 105 (22) (1996) 9982–9985.

[39] T. Zecho, A. Güttler, X. Sha, B. Jackson, J. Kuppers, J. Chem. Phys. 117 (2002) 8486.

[40] C. Leforestier, R.H. Bisseling, C. Cerjan, M.D. Feit, R. Friesner, A. Guldberg, A. Hammerich, G. Jolicard, W. Karrlein, H.-D. Meyer, N. Lipkin, O. Roncero, R. Kosloff, J. Comput. Phys. 94 (1991) 59.

[41] D. Kosloff, R. Kosloff, J. Comput. Phys. 52 (1983) 35.

[42] G.D. Billing, Dynamics of Molecule Surface Interactions, Wiley, New York, 2000.

[43] H. González-Herrero, E.C.-d. Río, P. Mallet, J.-Y. Veuillen, J.J. Palacios, J.M. Gómez-Rodríguez, I. Brihuega, F. Ynduráin, 2D Mater. 6 (2) (2019) 021004.