H₂ Thermal Desorption Spectra on Pt(111): A Density Functional Theory and Kinetic Monte Carlo Simulation Study

Abstract

Theoretical investigation of the static and kinetic behaviors of H and H₂ on metal surface plays a key role in the development of hydrogenation catalysts and new materials with high H₂ storage capacity. Based on the density functional theory (DFT) calculation of H and H₂ adsorption on Pt(111), H(a) adatom strongly interacts with surface Pt; while H₂ weakly adsorbs on Pt(111). H(a) adatoms stably occupy the face-centered cubic sites on Pt(111) which agrees with the experimental LERS observations. By using kinetic Monte Carlo (kMC) simulation, the qualitative effects of the kinetic parameters on the H₂ TDS spectra indicate that the H₂ desorption peaks shift to the low temperature with increasing pre-exponential factor and decreasing desorption barrier. Simultaneously, the desorption peaks shift downwards and broaden to two peaks with the increase of the lateral interaction energy among H(a) adatoms. Using the kMC simulation based on DFT calculation, the predicted H₂ TDS spectra are well consistent with the experimental ones. It unanimously proves that the two peaks of TDS spectra are derived from the lateral interactions among H(a). This work provides the intrinsic kinetics of H(a) and H₂ on Pt(111) at an atomic level, and gives insight into the development of hydrogenation catalysts.

1. Introduction

The adsorption and reaction behaviors of H and H₂ on metals are of great significance due to the wide application of hydrogenation reactions in heterogeneous catalysis [1,2,3,4,5,6,7]. It also provides basic information for the development of new materials with high H₂ storage capacity [8,9,10,11]. As a typical modeling system, the interactions between hydrogen and Pt(111) have attracted much interest and have been widely studied in experiment and theory [12,13,14,15,16,17,18,19,20].

Using high-resolution electron energy loss spectroscopy (HREELS), low-energy electron diffraction (LEED) and thermal desorption spectra (TDS), some consensus issues have been reached. H₂ dissociatively adsorbs on Pt(111) at low temperature; the maximum coverage of hydrogen (θ_H) is up to one monolayer (ML) below 100 K [21]; θ_H gradually decreases with the increase of temperature [22,23]; H adsorbs at the face-centered cubic (fcc) site verified by using low-energy recoil scattering (LERS) spectra [24,25]; diffusion activation energy of H on Pt(111) is much lower than that of H₂ desorption. These experimental conclusions were also supported by the related theoretical works [17,18,19,20,21,26,27,28]. Gudmundsdóttir et al. [17] investigated the H adsorption behaviors on Pt(110)-(1×2) by DFT calculations and found that the binding energy and activation energy for desorption are strongly dependent on hydrogen coverage. Thermodynamic description of hydrogen electroadsorption on Pt(111) was calculated by Hanh et al. [18] using first principles. Combining with a Monte Carlo simulation, the predicted cyclic voltammetry data is consistent with the experimental one. Hydrogen adsorption on Pt(100), (110) and (111) at high coverage was reported by Shi et al. [19] using thermodynamics based on DFT calculations. Their results indicated that hydrogen adsorption might modify the morphology of Pt catalysts. Karlberg et al. [20] proposed a cyclic voltammogram for H on Pt(111) and (100) from first principles. Combining experimental study of the H-D exchange reaction and first principles calculations, Xu [21] found that the repulsive interactions between H(a) increase with H(a) coverage, which leads to a reduction in both the adsorption energy and the desorption activation energy. Olsen et al. [26] found that a DFT/GGA approach with scalar relativistic effects is capable of describing the interaction between H(a) and Pt(111). The theoretical results by Watson et al. [27] indicated that H(a) diffusion barrier on Pt(111) is considerably smaller than that on Ni(111) and Pd(111) surfaces. By using the extended Hückel molecular orbital theory, Papoian [28] also found that H(a) moves freely on Pt(111).

Despite of the extensively studies of H/Pt(111) system, there are still some issues in debate. Christmann et al. [23] and Xu et al. [21] both reported that two H₂ desorption peaks at the low and high temperatures were observed in the TDS spectra on Pt(111). Since there is only one stable configuration for H adsorbed at the fcc sites on Pt(111) which has been proved by both experimental and theoretical results, what is the main reason leading to the two peaks? Which factors do determine the peak position and area of H₂ TDS spectra? Are the factors of H₂ desorption barrier (E_a), collision frequency among H(a) (pre-exponential factor, v), or the lateral interaction energy among H(a) on Pt(111) (E_int), etc.? What are the qualitative and quantitative effects for E_a, v and E_int on the H₂ TDS spectra? These important questions have still not been reported in the literature.

With the rapid development of multi-discipline cooperation in the study of the chemical reaction mechanism, kinetic Monte Carlo (kMC) simulation based on density functional theory (DFT) calculation has been applied in exploring the catalytic reaction mechanism [10,29,30,31]. The reaction thermodynamics and the potential energy surfaces of the possible reaction pathways on the modeling catalysts can be calculated by DFT theory [32,33,34,35,36]; while the kinetic reaction information on the specified surface can be obtained by the kMC simulation [30,37,38,39,40,41]. Hence, kMC simulation based on DFT calculation might provide the important information for the rational design of catalysts [42,43].

In this work, using DFT calculation, the adsorption, diffusion and desorption behaviors of H(a) and H₂ on Pt(111), as well as their electronic interactions are obtained. Then, the effects of H₂ desorption barrier, collision frequency and the lateral interactions among H(a) on the H₂ TDS spectra, are systematically investigated by kMC simulation. Finally, we combine the DFT calculation and kMC simulation to explore the intrinsic kinetics of the H₂ TDS spectra on Pt(111). This work provides important information regarding the interactions between H and metal materials, and sheds insight into the development of hydrogen storage material and hydrogenation catalysts.

2. Results and Discussion

2.1. DFT Study of H₂ and H on Pt(111)

We systematically investigated the adsorption behaviors of H and H₂ on Pt(111) by density functional theory calculation. The optimized structures are shown in Figure 1 and the corresponding parameters are listed in

Table 1. Parameters of H(a) and H₂ on Pt(111).

	Eads/kcal/mol	ZPVE Corrected Eads/kcal/mol	dPt-H/Å	dH-H/Å	Bader Charge/e
					H	Pt
H-fcc	−8.86	−9.12	1.874		−0.100	−0.027
			1.872			0.059
			1.875			−0.019
H-hcp	−7.47	−7.89	1.879		−0.057	−0.062
			1.866			0.036
			1.878			−0.004
H-top	−8.35	−7.35	1.572		−0.035	0.029
H2-fcc	−0.09	0.13	3.789	0.748
			3.712		0.025
			3.536		−0.032

. It can be learned from Table 1 that the adsorption energy of H(a) at the top site of surface Pt is −8.35 kcal/mol with the Pt-H bond length of 1.572 Å. The adsorption energies of H(a) at the fcc and hcp sites are −8.86 kcal/mol and −7.47 kcal/mol, respectively. The corresponding bond lengths of Pt-H are 1.866~1.879 Å. These results are well consistent with those reported by Graeme et al. [27] and Källén et al. [44].

Since the adsorption energy difference between H-top and H-fcc configurations on Pt(111) is just 0.51 kal/mol, the zero-point vibrational energies (ZPVE) are further calculated to determine which one is the most stable structure. The adsorption energies at the top and fcc sites are −7.35 and −9.12 kcal/mol, respectively. Therefore, H(a) at the fcc sites is the most stable configurations for H(a) on Pt(111), which is well consistent with the experimental results [24,25]. However, the adsorption energy of H₂ on Pt(111) with and without the ZPVE energy is 0.13 and −0.09 kcal/mol, respectively. The H-H bond length is 0.748 Å, which is almost the same as that of free H₂ molecule. Meanwhile, the distances of Pt-H are larger than 3.5 Å. This indicates the very weak interaction between H₂ and Pt surface. This will be further proved by the following analysis of the electronic properties of H(a) and H₂ on Pt(111).

To explore the electronic interactions of H(a) and H₂ with Pt(111), charge density differences (CDD) are calculated and shown in Figure 2. It is defined as: Δρ = ρ_{adsorbate/slab} − ρ_slab − ρ_adsorbate, where ρ_{adsorbate/slab} is the charge density for the adsorbed system; ρ_slab and ρ_adsorbate are the charge densities for the non-interacting slab substrate and adsorbate, respectively. The yellow and grey–blue parts represent the accumulation and depletion electrons, respectively.

As shown in Figure 2, H(a) atom gains electrons from the d orbital of surface Pt atoms; meanwhile the strong Pt-H bonds form due to the electron share between Pt and H(a). However, the electron transfer between H₂ and Pt(111) is almost negligible, which is consistent with the very weak interactions with the adsorption energy of −0.09 kcal/mol. This conclusion is also in agreement with the result of Bader charge as listed in Table 1. That is, H(a) gains electrons with negative charge; Pt loses electrons with partially positive charge.

As shown in the partial density of states (PDOS) of H(a) in Figure 3a, the s peaks of H(a) broaden and shift to −8~−6 eV for H-fcc and H-hcp configurations and −6~−4 eV for H-top in comparison with the free radical H atom. Meanwhile, there are the corresponding electronic peaks of Pt at the same energy position, as shown in the green peaks of Figure 3b. This indicates the strong interactions between H(a) and Pt(111). However, the PDOS peak shape of H₂ on Pt(111) remains unchanged, and the corresponding peaks of Pt are almost the same as those of the non-interacting substrate. Therefore, the very weak interactions between H₂ and Pt(111) occur.

2.2. kMC Simulation of H₂ TDS on Pt(111)

Before the kMC simulation of H₂ TDS on Pt(111) based on DFT calculation is performed, it is necessary to give a qualitative analysis of the effects of the main kinetic parameters on H₂ TDS on PT(111). This is very important for exploring the key parameters in kMC simulation.

A surface of 128 × 128 Pt atoms is adopted with the different H(a) coverages, in which H(a) atoms are randomly populated on Pt(111). The pre-covered H(a) models with different coverages mimic the experimental ones, in which different doses of gas atomic hydrogen are exposed to Pt(111) [21]. The samples with different H(a) coverages are heated from 100 K to 600 K with a heating rate of 10 K/s, which also agrees with the ultrahigh vacuum experimental conditions [21]. The H₂ TDS spectra include the following processes.

(1)

H adatom diffusion: H(*) + * → * + H(*)

(2)

H₂ desorption: 2H(*) → H₂(g) + 2*

(3)

H₂ dissociation: H₂(g) + 2* → 2H(*)

The asterisk * with and without species denotes an adsorbed species and empty surface site, respectively. Since the Leybold stainless-steel UHV chamber is continuously vacuuming with the total pressure of about 10⁻⁸ Pa [21], the desorbed H₂ molecules will be quickly pumped out of the UHV chamber. Therefore, the dissociation of H₂ is not involved in the kMC simulation.

To explore the intrinsic kinetics of H₂ TDS spectra on Pt(111), we investigate the effects of the different kinetic parameters on the spectra, as shown in Figure 4, Figure 5 and Figure 6.

2.2.1. Effects of the Pre-Exponential Factor of H₂ Desorption

As shown in Figure 4, one symmetric desorption peak is observed in the H₂ TDS spectra at the low (0.12 ML), medium (0.6 ML) and high (1.0 ML) coverages, respectively. The peaks shift to the low temperature with the increase of pre-exponential factor (v). Since the pre-exponential factor represents the collision frequency of reactants upon reaction, the higher v leads to the more possibility of H(a) desorption at the same temperature. Therefore, the peak of H₂ TDS shifts downwards with the pre-exponential factor. Meanwhile, the peaks become higher and narrower than those of v = 10¹² s⁻¹.

2.2.2. Effects of the Desorption Barrier of H₂ Desorption

Similar results are observed for the H₂ TDS spectra on Pt(111) with the decrease of activation energy for H(a) desorption (E_a). That is, only one desorption peak is observed in Figure 5, and it shifts downwards with the decrease of E_a. Because the H(a) desorption rates increase with the low desorption barrier at the same temperature, the peak of H₂ TDS spectra shifts downwards with the decrease of E_a. Similarly, the corresponding desorption peaks also become higher and narrower than those of E_a = 24 kcal/mol.

2.2.3. Effects of the Lateral Interactions among H(a)

To systematically investigate the effects of the lateral interactions among the nearest-neighbor H(a) on the H₂ desorption behaviors, the H₂ TDS spectra on Pt(111) with the different lateral interaction energies at the low, medium and high coverages are simulated and shown in Figure 6.

At the low coverage (0.12 ML), the H₂ desorption peaks at 390 K almost remain unchanged with increasing the lateral interaction energy (E_int). This might be attributed to the large distance among H(a) species at the low coverage, which leads to the invalid lateral interaction effects on the TDS spectra. If the H(a) coverage increases to 0.6 ML, the main desorption peaks shift to 350, 320 and 290 K corresponding to the E_int of 0.2, 0.4 and 0.6 kcal/mol, respectively. Meanwhile, two peaks appear when E_int is larger than 0.2 kcal/mol. When the coverage is up to 1.0 ML, the highest desorption peaks shift to 330, 290 and 260 K with E_int, respectively. This is attributed to the strong effects of the lateral interactions, which leads to the decrease of the desorption barrier of H(a) species.

2.3. kMC Simulation of H₂ TDS on Pt(111) Based on DFT Calculation

2.3.1. Parameters Calculated by DFT Method

According to the qualitative analysis of the effects of main kinetic parameters on H₂ TDS spectra, it indicates that the activation energies of H₂ desorption, pre-exponential factors and the lateral interaction energy of H(a) play the important roles in the kMC simulation of H₂ TDS on Pt(111). Hence, the key parameters should be calculated by the DFT method before the kMC simulation based on first principles.

The predicted barriers of H₂ desorption and H(a) diffusion on Pt(111) are shown in Figure 7. It can be drawn from Figure 7 that the desorption barrier of H₂ on Pt(111) is 21.29 kcal/mol. For the transition state structure, two H(a) atoms move to the top site of surface Pt atom with the Pt-H bond lengths of 1.691 and 1.754 Å, respectively. The distance between the two H(a) is 0.967 Å. If they are further close to each other, H₂ forms and desorbs from Pt(111) surface. The image frequency of TS structure is 377 cm⁻¹, which its vibration direction points to the dissociated H(a) and H₂, respectively. The diffusion barrier of H(a) at the fcc site to hcp site is 1.71 kcal/mol via a bridge site with the H-Pt bond lengths of 1.832 Å.

Based on the harmonic transition-state theory, the pre-exponential factors (ν_i) of each reaction step are calculated using the following equation: $v_i=\frac{{\displaystyle {\prod }_1^{3N}{f}_i^{\mathrm{IS}}}}{{\displaystyle {\prod }_1^{3N-1}{f}_i^{\mathrm{TS}}}}$, where $f_i^{\mathrm{IS}}$ is the vibrational frequency at the initial state, and $f_i^{\mathrm{TS}}$ is the vibrational frequency at the transition state (excluding the imaginary one). For the reverse pathway, the vibrational frequencies at the final state ($f_i^{\mathrm{FS}}$) are used. Frequencies used to calculate the pre-exponential factors are listed in Tables S2–S4 of the supplementary materials.

It is reported that the lateral interactions among H(a) are mainly derived from the short-range pair interaction by the nearest neighbor (NN) adsorbates through surface, and the interactions among the next nearest neighbor (NNN) H(a) adatoms are generally negligible [21,45]. Therefore, the lateral interactions among the NN H(a) are calculated by DFT theory in this work. Since the lateral interactions among H(a) atoms are relatively weak, four models shown in Figure 8 are designed to obtain the accurate value. The lateral interaction energy is calculated by the following equation: $E_{\mathrm{int}}=\frac{1}{m}\left(\frac{1}{n}{E}_{\mathrm{nHads}}-E_{1\mathrm{Hads}}\right)$, where n is the number of H(a); m is the number of pairwise interactions per H(a) atom. The corresponding parameters are listed in

Table 2. Parameters of lateral interaction energy among H(a) adatoms.

	Eads/kcal/mol	Eav/kcal/mol	Eint/kcal/mol	m
1H-fcc	−8.86	−8.86	-	0
2H-fcc	−16.93	−8.46	0.39	1
3H-fcc	−24.49	−8.16	0.35	2
4H-fcc	−31.41	−7.86	0.40	2.5
9H-fcc	−61.55	−6.85	0.34	6

Note: m is the number of pairwise interactions per H(a) atom.

.

It can be learned from Table 2 that the average adsorption energies of H(a) decrease with coverage, which is derived from the lateral interactions among H(a). The number of pairwise interactions per H(a) atom for the nH-fcc (n = 2, 3, 4 and 9) on Pt(111) are 1, 2, 2.5 and 6, which corresponds to the lateral interaction energies are 0.39, 0.35, 0.40, 0.35 kcal/mol, respectively. Therefore, the average value of 0.37 kcal/mol is adopted in the following kMC simulation based on DFT calculation.

2.3.2. kMC Simulation Based on DFT Calculation

For the kMC simulation based on DFT calculation, all the parameters should be used by the theoretical predictions. However, the diffusion barrier of H(a) on Pt(111) is just 1.71 kcal/mol, indicating that most of computer time will be consumed in the calculation of rate constant by using the Arrhenius equation due to the very low diffusion barrier. To improve the efficiency of kMC simulation, the effects of H(a) diffusion rates are tested and shown in Figure S1 in the supplementary material. Based on the tests, the diffusion rate constant has no effects on the kMC simulation results if it is larger than 10 s⁻¹. In this work, the rate constant of 1000 s⁻¹ is adopted.

The H₂ TDS spectra calculated by the kMC simulation based on DFT method at different coverages are shown in Figure 9. It can be learned from Figure 9 that the desorption peak at 0.05 ML is at 345 K, which is well consistent with the experimental value of 343 K [21]. It shifts to the low temperature with the increase of H(a) coverage. When the coverage is larger than 0.6 ML, two peaks are observed. Furthermore, the desorption peak at the low temperature shifts to 250 K when the coverage is up to 1.0 ML. This is also consistent with the experimental peak temperature at 235 K in TDS spectra [21]. The H₂ TDS spectra derived from the kMC simulation based on DFT method unanimously prove that the two desorption peaks are derived from the lateral interactions among H(a) species. The simulation results are well consistent with those in experiment, indicating the method of kMC based on DFT calculation is feasible and reliable.

3. Model and Method

3.1. DFT Calculation Details

Spin-polarization calculations were performed using Vienna ab initio simulation package (VASP) code based on a first-principles density functional theory [46,47]. The electron exchange and correlation were treated within the generalized gradient approximation (GGA) of the revised Perdew Burke Ernzerhof (RPBE) functional [48]. The Kohn–Sham one-electron wave functions were expanded in a planewave basis with a cutoff energy of 400 eV. The Brillouin-zone integration was approximated by a sum over special k-points using the Monkhorst-Pack (MP) [49] grids and Methfessel–Paxton smearing method with a sigma value of 0.1 eV. The electrons in Pt 6s, 5d and H 1s orbitals were treated as valence electrons. The ionic cores were represented by the projector augmented wave (PAW) potentials [50]. The Kohn–Sham equations were solved self-consistently. The self-consistent field tolerance for the energy calculation was set as 1.0 × 10⁻⁴ eV; the maximum Hellmann–Feynman force tolerance for the structure optimization was set as 0.02 eV/Å.

The calculated bulk lattice constant (3.997 Å) of Pt(111) with a MP grid of 13 × 13 × 13 agrees well with the experimental value (3.924 Å) and the values reported in the literature [27,28]. A Pt(111) slab model with five layers was chosen for the theoretical calculations. The two bottom layers were fixed during geometric optimization. The vacuum space was set as 12 Å along the (111) direction to minimize the interaction between slab surfaces. The p(3×3) (R = 60°) lateral cell was involved in the present work with the k-points grid of 3 × 3 × 1. By testing the convergence of k-points, cutoff energy and size of Pt(111) slab, these settings ensure the required precision.

The total adsorption energy of nH(a) on Pt(111) was defined as E_ads=E_nH/slab − (E_slab + 0.5nE_H2), in which E_(nH/slab), E_(slab), and E_(H2) represent the total energies of a slab model with adsorbed hydrogen atoms, the slab, and free hydrogen molecule, respectively. The average adsorption energy per H(a) on Pt(111) was calculated by E_av = E_ads/n. The adsorption energy of H₂ was defined as E_ads = E_H2/slab − (E_slab + E_H2). The vibrational properties were calculated on the basis of the harmonic vibration of the adsorbed H(a) and H₂(a) on Pt(111). The zero-point vibrational energy (ZPVE) is calculated by the equation of $E_{\mathrm{ZPVE}}=\frac{1}{2}{\displaystyle \sum _{i=1}^n{f}_i}$, where f_i is the energy of ith vibrational mode, and its unit should be converted from cm⁻¹ to kcal/mol or eV. The frequencies used in ZPVE corrections were listed in Table S1. Transition states (TS) of H(a) desorption and diffusion were calculated with the climbing image nudged elastic band (CINEB) method [51]. The activation barrier was defined as the total energy difference between the transition state and the corresponding stable structures of reactants or products.

Since it was reported that the van der Waals interactions between adsorbates and surfaces do influence the adsorption and reaction behaviors of reactants and intermediates [52,53], we also compared the adsorption energies of H(a) and H₂ on Pt(111), as well as the barriers of H(a) diffusion and H₂ desorption using DFT and DFT-D methods with RPBE and PBE functionals, as listed in Table S5. The theoretical results by the two methods with PBE functional are almost the same. This might be due to two facts. One is that the volume of hydrogen atom is very small, which leads to the negligible interactions between nH-fcc through space. The other is that the lateral interactions of nH-fcc are mainly derived from the repulsive ring among H(a) through Pt(111) surface, which the corresponding electronic repulsive interactions can be well calculated by DFT method. Therefore, the results calculated by DFT method are reliable.

3.2. kMC Simulation

The kMC simulation was carried out on a periodic p(128×128) (R = 60°) superlattice of Pt(111) surface, which represents the single crystal surface of Pt(111) in the ultrahigh vacuum chamber. According to the experimental results, H atoms only adsorb at the face-centered cubic (fcc) site of Pt(111); while H₂ very weakly physisorbs on Pt(111). Hence, the initial model of H(a) on Pt(111) with a designated coverage was produced by H random distribution at the fcc sites on Pt(111). The temperature-programmed desorption rate of H at Pt(111) was set as 10 K/s. The rate constants of the elementary surface step follow the Arrhenius equation: r_D = vexp(−E_a/kT), where v is the pre-exponential factor; T is the temperature, and E_a is the activation barrier for the elementary processes [40,54]. The quantitative effects of the lateral interactions on the H₂ desorption barrier is described by the effective activation energy defined as: E_eff = E_a − ∑xE_int, where E_a is the activation energy of H(a) desorption at zero coverage; x is the number of nearest-neighbor H(a); E_int is the lateral interaction energy [31]. According to the definition, the effects of the lateral interactions on the transition state originate from an interpolation between the initial and final states due to the repulsion among the nearest neighbor H(a). The probability of a process to be chosen and implemented is directly proportional to the rate calculated by the master equation: [40] $\frac{\mathrm{d}{P}_{\mathsf{\alpha }}}{\mathrm{d}t}={\displaystyle \sum _{\mathsf{\beta }}\left(k_{\mathsf{\alpha }\mathsf{\beta }}{P}_{\mathsf{\beta }}-k_{\mathsf{\beta }\mathsf{\alpha }}{P}_{\mathsf{\alpha }}\right)}$, in which α and β are the adsorbates at the different sites; P_α and P_β are their probabilities of forming α and β; k_αβ is the rate constant of α transforming into β. The master equation was solved by the discrete event simulation algorithm. There are many kMC simulation programs, such as CARLOS [55], Multi-scale KMC [56] and kmos [57], etc. The program herein was written in C language [58]. The first reaction method with the time-dependent rate constants was adopted in the work. A tentative time was calculated for every possible process [42]. Consequently, the numbers of steps of the temperature programmed desorption simulation could be obtained.

4. Conclusions

According to the DFT calculation of H and H₂ adsorption on Pt(111), H(a) at the fcc site is the most stable structure. Strong interactions between H(a) and surface Pt occur; while H₂ weakly physisorbs on Pt(111). According to the qualitative effects of the kinetic parameters on the H₂ TDS spectra, it is found that the desorption peaks shift to the low temperature with increasing v and decreasing E_a of H₂ desorption. Simultaneously, the desorption peaks become higher and narrower. With the increase of E_int, the desorption peaks shift downwards and broaden to two peaks. By the kMC simulation based on DFT method, the predicted H₂ TDS spectra are well consistent with the experimental ones. This work provides the intrinsic kinetics of H(a) and H₂ on Pt(111) at the atomic level, and gives insight into the development of hydrogen storage material and hydrogenation catalysts.