Isotope Effect on Hydrogen Desorption from a Graphene Flake by Ab Initio Electron – Ion Dynamics Simulation

Thermal and femtosecond-laser-assisted desorption of hydrogen (H) and deuterium (D) from a graphene flake were investigated by performing ab initio electron–ion dynamics simulations. In thermal desorption, the isotope effect appears only in the desorption velocity of H and D atoms. In contrast, the isotope effect is pronounced in laser-assisted desorption. H atoms are desorbed but D atoms remain adsorbed when irradiated by a laser pulse with a certain power and frequency. The reason is that the laser pulse excites electronic states that critically modify the forces on H (D) atoms and carbon (C) atoms to which H (D) atoms are bound. The present simulation predicts that it is possible to precisely control the desorbing species by tuning the laser parameters. [DOI: 10.1380/ejssnt.2011.1]


I. INTRODUCTION
Since the isolation of graphene [1,2], a single atomic layer of hexagonal carbon, research on graphene has been increasing enormously because of its remarkable physical properties and a stable atomic structure that makes it chemically inert.However, its unique electronic structure that has a gapless energy band prevents graphene from being applied to nanoscale electronic devices with on and off states.Thus, there have been various efforts to open the electronic band gap [3,4].A recent experiment on hydrogenated quasi-free-standing graphene achieved a tunable band gap ranging from 1 to 3.5 eV by varying the hydrogen (H) coverage [5].Therefore, adsorption/desorption of H on graphene are potentially key techniques for realizing graphene-based electronic devices.
First-principles electronic-structure calculations clarified the energetics of H adsorption on graphene [6][7][8][9].The dynamic properties of H and deuterium (D) on graphite surfaces have been studied experimentally for thermal desorption [10] and femtosecond-laser-assisted desorption [11].Interestingly, both experiments clearly revealed an isotope effect in the desorption spectra.An isotope effect has also been observed in time-of-flight spectra of photodesorption of H physisorbed on Ag and it was explained by coherent desorption induced by electronic transitions model [12].Although theoretical analysis based on a model Hamiltonian is useful for interpreting spectra [11,12], ab initio dynamics simulations are required to determine the real-time H (D) desorption mechanism.In particular, since laser-assisted desorption is triggered by ultrafast electronic excitations, the correlated electronion dynamics has to be treated from first principles [13].
The objectives of the present study are to determine the real-time desorption dynamics of H and D in thermal and laser-assisted desorption of graphene and to interpret the isotope effect on laser desorption using timedependent density functional theory (DFT) in combination with molecular dynamics (TDDFT-MD) simulations [13,14].We chose pyrene (C 16 H 10 ) as a simple model for graphene to reduce the computational time of the present simulation.The dynamics of the present model is thought to reveal the intrinsic behavior occurring locally on graphene.The remainder of this paper is organized as follows.In §II, we describe the techniques of DFT and TDDFT-MD and a model for H (D) adsorbed pyrene C 16 H 10 .In §III, we present the results for potential energy surfaces (PESs) of the H atom and the time evolution of H (D) and C atoms for thermal desorption and laser-assisted desorption.We discuss the critical role of C atoms to which H (D) atoms are bound in desorption.Finally, §IV gives concluding remarks.

A. Methods
We briefly describe the DFT, DFT-MD, and TDDFT-MD schemes employed in the present calculations.We first used conventional DFT in real space to calculate the potential energy curves for an H atom binding to a pyrene molecule (C 16 H 10 ) in the ground state.The computational details are as follows.The Kohn-Sham equation is solved self-consistently using the higher-order finitedifference method [15] with a nine-point formula.The mesh spacing is 0.159 Å.The calculation box has dimensions 8.47 × 9.74 × 11.0 Å3 .The interactions between the ionic cores and the electrons are described by a normconserving pseudopotential [16] formulated by Troullier and Martins [17] based on the Kleinman-Bylander separable form [18,19] in real space.The exchange-correlation energy is treated within the local-spin-density approximation (LSDA) using the Ceperley-Alder [20] functional with the Perdew-Zunger [21] parameterization.The spinpolarization effect is found to be important for correctly describing the long-range correlation energy [6].
We then performed DFT-MD calculations to simulate the atom dynamics for thermal desorption of H and D atoms from C 16 H 10 .The initial kinetic energy of 1.5 eV is imparted to H and D atoms and the time evolution of the positions of all atoms are determined by solving the classical Newton's equation of motion with the Hellmann-Feynman (HF) forces on the atoms.The HF forces are calculated under the condition of the electronic ground states of the whole system at each MD time step.Namely, the present dynamics for thermal desorption is on the Born-Oppenheimer surface.The time step of the DFT-MD calculation is 4.84 ×10 −2 fs.
We finally performed TDDFT-MD simulations [14,22] to determine the time evolutions of the electronic wave functions and the atomic positions after H (D) adsorbed C 16 H 10 was irradiated by a femtosecond laser pulse.The time evolution of the electronic wave functions is determined by solving the time-dependent Kohn-Sham equations with the adiabatic LSDA for the exchangecorrelation functional.The time-development operator is approximated by a fourth-order Taylor expansion [23].The atomic positions are determined by solving the classical Newton's equation of motion with similar HF forces on the atoms as those in the thermal desorption dynamics described above.Since the TDDFT-MD scheme may give rise to transitions to some averaged states of the adiabatic potential energy surfaces, the atoms follow Ehrenfest dynamics [24].The laser pulse is approximated by a spatially uniform electric field that has a Gaussian temporal profile (see Fig. 4(a)) [13].The time step employed in the time-dependent calculations is 4.84 × 10 −4 fs.Refer to Ref. [13] for other computational details.

B. Models
The present system, a pyrene molecule (C 16 H 10 ) with an H (D) adatom (see Fig. 1) was selected as a minimal model for investigating the desorption dynamics of single atoms from graphene.The 10 hydrogen atoms serve only to passivate the terminal dangling bonds of the C 16 cluster.Although this cluster appears to be too small to reproduce the binding properties of an H (D) atom to a large graphene sample, the calculation results agree well with those of other numerical studies that employ a larger cluster (C 24 H 12 ) [6] or a system with periodic boundary conditions and a (2 × 2) unit cell [8].The C-C and C-H bond lengths of an isolated C 16 H 10 cluster are 1.40 and 1.09 Å, respectively.Figures 1(a

A. Potential energy surfaces for H on C16H10
PESs similar to those of the present results have been reported by other groups [6,8], although they have slightly different C-atomic planes from ours (see § II B).Here, we focus on only the stable position of C to which H (D) is bound because the C atom plays an important role in H and D desorption, especially for laser-induced desorption, which is discussed in § III C. In contrast to H on metal surfaces, H chemisorbed on graphene does not diffuse in the plane because of the high potential barriers parallel to the plane [8].We thus calculated the PESs for an H atom along the direction normal to the C 16 H 10 plane by assuming low H coverage and no recombinative desorption.
Figure 2(a) shows the potential energies of an H atom on C 16 H 10 as a function of Z H .The broken curve is for the PES when the C position is fixed, whereas the solid curve is for the PES obtained when the C position is relaxed during structure optimization.Figure 2(b) shows the C position with regard to Z H for two optimization conditions for the C atom (broken and solid curves).Interestingly, the C atom initially follows the H atom before suddenly detaching from the H atom and returning to its position in the C 16 H 10 plane when Z H exceeds ∼2 Å (see Fig. 2(b)).The discontinuous change in the C position reflects the kink in the PES (solid curve in Fig. 2(a)).Consequently, it is reasonable to expect that the H (D)desorption dynamics depends critically on the kink in the potential, i.e., on C detachment.This is why we present the C-position-dependent PESs before discussing the isotope dependence of the desorption dynamics.http://www.sssj.org/ejssnt(J-Stage: http://www.jstage.jst.go.jp/browse/ejssnt/) e-Journal of Surface Science and Nanotechnology

B. Thermal desorption dynamics
We present the numerical results for the time evolutions of H, D, and C atoms obtained by performing DFT-MD simulations.The following initial conditions were used.The kinetic energy of 1.5 eV is imparted only to H and D atoms adsorbed on top of the C atom puckered from pyrene in the direction normal to the plane.We chose the simplest initial conditions for thermal desorption to demonstrate the different motions of H and D atoms and the key role of the C atom in real time.All atoms move in accordance with the forces acting on them.
Figure 3 shows the time evolutions of the positions of the H, D, and C atoms.The blue and red curves in Fig. 3(a) represent Z H and Z D , respectively.In addition, the blue (red) curves in Fig. 3 However, this isotope effect on the desorption dynamics obtained using the present simulation does not explain the different activation energies for H and D desorptions [10] that have been observed experimentally.The experimental activation energies for H and D desorptions are 0.6 and 0.95 eV, respectively.These different activation energies cannot be attributed to the different zero-point vibration energies (ZVEs) of H and D because the difference in the ZVEs is evaluated to be only about 0.03 eV from the potential curve in Fig. 2(a), which is one-tenth the difference in the activation energies observed in the experiment.
We consider that our simulation does not predict the observed difference in the activation energies because of the artificial initial conditions used for the cluster atoms.Namely, only desorbing H or D atoms are assumed to have initial kinetic energies in the simulation, instead of all the atoms in the pyrene molecule.Therefore, to allow direct comparison with experimental observations, the initial conditions used when simulating thermal desorption should consider the thermal motions of all atoms at the experimental temperature.

C. Laser desorption dynamics
In this section, we present the results for the desorption dynamics of H and D atoms induced by irradiation by an intense ultrafast laser pulse.We used the TDDFT-MD simulation to account for the effect of forces induced by electronic excitation on desorbing atoms.We applied a laser pulse with an electric field strength of 6.7 V/ Å, a frequency of 2.7 eV, and a pulse width of 30 fs (see Fig. 4(a)).The laser pulse was polarized perpendicular to the cluster plane.
We found that H and D atoms exhibit substantially different time evolutions (see Fig. 4(b)) and that the C atom has different time evolutions for H and D desorption (see Fig. 4(c)); these different time evolutions are in marked contrast with those for thermal desorption (see Fig. 3).Both H and D atoms move toward the vacuum while the C atom moves toward the cluster plane up to 20 fs when the laser field reaches its maximum intensity.After 30 fs, when the laser field ceases, the C atom returns toward the vacuum and captures a D atom.Consequently, the D atom is pulled toward its original position in the cluster.However, the C atom cannot capture an H atom and it continues to move toward the cluster plane beyond t = 20 fs, because the H atom moves to the vacuum faster than the D atom.Thus, H is desorbed but D is not desorbed.The complex features in the time evolutions of these atoms are attributed to the time-varying direct force of the laser electric field acting on the atoms and to an indirect force caused by potential generated by laser-induced electronic excitations [13].
Therefore, the desorption dynamics tends to be irreg-ular and critically depends on the laser parameters.Actually, both H and D atoms were desorbed when a laser pulse with an field of 7 V/ Å was used, which is marginally higher than that used above (6.7 V/ Å).Further, we found that the motion of the C atom to which the H (D) atom is bound gives rise to unusual desorption dynamics.Here, we mention the laser power, 5×10 14 W/cm 2 in the present simulation.The value is higher than the typical one of femtosecond laser pulses, ∼ 10 13 W/cm 2 in experiments [25].Our high power laser, however, cannot fully excite the electronic states of the pyrene cluster, because the frequency of present laser does not corresponds to the one of the particular frequencies of the laser at which the photoabsorption [13] by the pyrene predominantly occurs.As a result, the total energy transfer from the applied laser to the pyrene is very much reduced to be ∼ 30 mJ/cm 2 , which is a typical order of magnitude in the experiments.This indicates that we could extract D atoms with the laser of much lower power by tuning the laser frequency to enhance the photoabsorption of the pyrene.Such detailed simulations will be carried out near future.The present TDDFT-MD simulations indicate that the isotope effect is pronounced for the desorption dynamics under laser irradiation.The origin of this enhanced isotope effect is the time-varying and excited-state forces induced by ultrafast laser fields that simultaneously give rise to transitions to some averaged states of the adiabatic PESs [13].In contrast, there is a rather weak isotope effect on thermal desorption dynamics (see Fig. 3) because H and D atoms experience the same adiabatic PES for thermal desorption.
The isotope effect on femtosecond-laser desorption of H and D atoms has been observed and analyzed using the Newns-Anderson model by Frigge et al. [11].It is not possible to directly compare their results with the present results because they employed different adsorption geometries for H and D atoms (ortho or para positions), different laser parameters, and a different methodology.However, the isotope characteristics obtained in the present study (see Fig. 4(b)) are consistent with the velocity distributions of desorbing H and D atoms obtained in their study.Specifically, they found that the total desorption yield is much higher for H atoms than for D atoms and that the velocity that gives the highest yield is larger for H atoms than for D atoms.

IV. CONCLUSION
We investigated thermal desorption and laser-assisted desorption dynamics of H and D atoms from H (D) adsorbed C 16 H 10 by performing DFT-MD and TDDFT-MD simulations.We found the isotope effect in both desorption processes.The isotope effect is not pronounced for thermal desorption.Specifically, with the exception of the desorption velocity (which is higher for H atoms than for D atoms), H and D atoms have similar dynamics.In contrast, the isotope effect is pronounced for laser desorption.H atoms are desorbed but D atoms remain adsorbed when irradiated by a laser pulse with a certain power and frequency.The forces generated by the time-varying elec-tric field and the excited electronic states cause H and D atoms to have very different time evolutions.Significantly, the motion of the C atom to which a H (D) atom is bound is found to be irregular and critically influences the dynamics of the adsorbates.The present simulation not only reveals the femtosecond dynamics for desorption in real time but it also implies that desorption can be controlled by tuning the laser parameters.

FIG. 1 :
FIG.1: H-adsorbed pyrene (C16H10) as a model cluster of an H-adsorbed graphene.C and H atoms are depicted by white and blue circles, respectively.(a) Top and (b) side views of the model cluster.ZH is the distance between adsorbed H and the surface plane and ZC is the distance between puckered C and the surface plane.The initial velocity for the thermal desorption dynamics ( §III B) and the direction of the electric field for laser desorption dynamics ( §III C) are respectively denoted by v and E.
) and (b) re-spectively show top and side views of the stable structure of H-adsorbed C 16 H 10 .The H sits stably on top of the C atom, pulling the C atom up toward the vacuum by a distance of 0.3 Å from the C 16 H 10 plane.Z H and Z C respectively refer to the distances of the H and C atoms from the C 16 H 10 plane.

FIG. 2 :FIG. 3 :
FIG.2: (a) Calculated potential energy curves for an H atom above a C atom in a pyrene cluster.The solid curve (closed circles) is obtained when the puckered C atom is allowed to fully relax for each H position.The broken curve (open circles) is obtained when all the surface atoms are fixed to their equilibrium positions.(b) ZC as a function of ZH.When ZH exceeds 2 Å, the C-H bond suddenly breaks and the C atom returns to its unperturbed surface position.The conditions of calculation for solid and broken curves are the same as those of (a).
(b) indicate the time evolutions of Z C when H (D) atoms are desorbed.The H and D atomic motions depicted in Fig. 3(a) differ only in terms of the desorption velocity.This is expected since the H and D atoms had the same kinetic energies.The different desorption velocities of the H and D atoms induce differ-ent C atom motions because of the different interaction times of H-C and D-C.The peak value of Z C is larger for D desorption than for H desorption in Fig. 3(b).This is probably because H and D desorptions give rise to different motions for the other atoms in the cluster; otherwise the peak position of Z C for D desorption should shift to the right and the Z C peak height should remain the same as that for H desorption.

FIG. 4 :
FIG. 4: Laser desorption dynamics.An H-adsorbed pyrene cluster is irradiated by a laser pulse that is polarized perpendicular to the cluster plane.(a) Profile of the electric field of the laser pulse.(b) Time evolution of ZH.The blue (red) curve represents the H (D) position.(c) Time evolution of ZC.The blue (red) curve represents the C position for H (D) desorption.