Understanding the Role of Rutile TiO2 Surface Orientation on Molecular Hydrogen Activation

Titanium oxide (TiO2) has been widely used in many fields, such as photocatalysis, photovoltaics, catalysis, and sensors, where its interaction with molecular H2 with TiO2 surface plays an important role. However, the activation of hydrogen over rutile TiO2 surfaces has not been systematically studied regarding the surface termination dependence. In this work, we use density functional theory (PBE+U) to identify the pathways for two processes: the heterolytic dissociation of H2 as a hydride–proton pair, and the subsequent H transfer from Ti to near O accompanied by reduction of the Ti sites. Four stoichiometric surface orientations were considered: (001), (100), (110), and (101). The lowest activation barriers are found for hydrogen dissociation on (001) and (110), with energies of 0.56 eV and 0.50 eV, respectively. The highest activation barriers are found on (100) and (101), with energies of 1.08 eV and 0.79 eV, respectively. For hydrogen transfer from Ti to near O, the activation barriers are higher (from 1.40 to 1.86 eV). Our results indicate that the dissociation step is kinetically more favorable than the H transfer process, although the latter is thermodynamically more favorable. We discuss the implications in the stability of the hydride–proton pair, and provide structures, electronic structure, vibrational analysis, and temperature effects to characterize the reactivity of the four TiO2 orientations.


Introduction
Titanium oxide (TiO 2 ) has been widely used in numerous fields, from everyday applications (paint, inks, toothpaste, makeup) to technological devices, such as dye-sensitized solar cells (DSSCs) [1,2], photoelectrochemical cells [3], photocatalysts [4], catalysis [5,6], sensors [7,8], biomedical treatments [9], lithium ion batteries [10], or photovoltaics [11,12]. The interaction of hydrogen with TiO 2 surfaces plays an important role in many reaction processes [13][14][15][16][17][18][19][20] and has been widely studied [21][22][23][24][25][26][27]. Despite the high interest generated by hydrogen-titania interfaces, the nature of the species involved is still poorly understood-protons are generally reported as being stable in hydrogenated rutile (110) [28], atomic surface hydrogen has been found to prevent electron-hole recombination on an Au-TiO 2 photocatalyst [14], and very recently hydride species have been characterized as being stable on its surface [29,30]. In this work, we investigate the role of the surface termination in the H 2 dissociation and migration on rutile surfaces. We focus on the characterization of the stability of surface Ti-H species formed by interaction with H 2 and their subsequent transfer to neighboring oxygen sites in order to provide a comprehensive picture of the adsorption, desorption, and diffusion mechanisms occurring at H-TiO 2 interfaces.

Materials and Methods
Density functional theory (DFT) calculations were performed using the Vienna ab initio simulation package (VASP) version 5.4.4 [68]. Projector-augmented wave (PAW) pseudopotential was used to describe the core electron representation with 1, 4, and 6 valence electrons for H, Ti, and O, respectively [69,70]. The generalized gradient approximation (GGA) approach was used for the exchange and correlation potential with the Perdew-Burke-Erzenhof (PBE) functional [71,72]. The GGA+U approach of Dudarev et al. was used to treat the 3d orbital electrons of Ti with the effective Hubbard on-site Coulomb interaction parameter (U' = U − J) [73]. We chose U' = 4 according to the proposed value from previous works [24,28,74], referred herein as U. A 400 eV cutoff energy for the plane-wave basis set was found to correctly treat the rutile surface [28]. The dissociation of hydrogen on rutile TiO 2 surfaces was investigated in the 1 × 1 unit cell for (001), (100), and (101) and in the 2 × 2 unit cell for the (110) surface. The open shell systems were treated with spin polarized calculations. The energy convergence was set to 3.0 × 10 −2 eV for the ionic loop and 1.0 × 10 −4 eV for the electronic loop. The slab models were cut from the optimized structure of bulk rutile (Figure 1). A vacuum layer of 20 Å was employed. The slab thickness used is given in Table 1. The lower-half layers of the slab were kept frozen and the upper-half layers were allowed to relax. We used the Monkhorst−Pack scheme to sample the Brillouin zone, and the distance between each k-point was 0.033 Å −1 [18,35]. The constrained minimization and climbing-image nudged elastic band (CI-NEB) methods were used to locate transition states (TS) [75,76]. In this work, the minimum energy pathway for each elementary reaction was discretized by a total of four images between the initial and final states. The imaginary frequency of every transition state was checked to connect initial and final states. The zero point energy (ZPE) vibration energy was calculated from vibrational frequencies as one-half of the sum of real-valued harmonic vibrational frequencies [77].  We also consider the effect of temperature by calculating the Gibbs free energy at room temperature (298 K); in the solid system, the pressure volume term pV can be ignored, thus:

G(T) = H − TS = U + pV -TS ≈ U(T) -TS(T)
It is reasonable to only consider the vibrational contributions, therefore:

S(T) = Svib(T)
For vibrational spectra, the density-functional perturbation theory (DFPT) linear response approach was used [78,79]. The matrix of Born effective charges (BEC) is obtained and indicates the change of involved atom's polarizabilities. The infrared intensity can be described as in the following formula containing Born effective charges and the eigenvectors ( ) e s β υ : where and β are Cartesian polarization, ( ) e s β υ indicates the normalized vibrational eigenvector, and Z αβ * indicates the effective charge tensor. To assess how the frequencies obtained depend on the computational setting, the performance of four different density functionals (PBE, Local Density Approximation (LDA), Perdew-Wang (PW91), Perdew-Burke-Ernzerhof revised for solids (PBESOL)), cut-off (300, 400, 500, 600, and 700 eV), choice of U (3, 4, 5, 6, and 7 eV), and the inclusion of dipole corrections were tested (see Supplementary Tables S1-S4 and Figures S7-S10). Although the numerical values are affected by computational settings, the trends between the different orientations are maintained.
Dispersion effects were evaluated for the heterolytic path for the intermediates and TS1 (the latter as a single-point calculation) by means of dispersion corrections (Grimme D3) (zero) [80] and the results are displayed in Supplementary Table S7 and Figure S12. Dispersion corrections were   (4) Ti (5) Ti (5,6) Ti (4,5) Parameter: a, b in Å a = 4.661 a = 4.661 a = 6.018 We also consider the effect of temperature by calculating the Gibbs free energy at room temperature (298 K); in the solid system, the pressure volume term pV can be ignored, thus: It is reasonable to only consider the vibrational contributions, therefore: For vibrational spectra, the density-functional perturbation theory (DFPT) linear response approach was used [78,79]. The matrix of Born effective charges (BEC) is obtained and indicates the change of involved atom's polarizabilities. The infrared intensity can be described as in the following formula containing Born effective charges and the eigenvectors e β (s|υ ): where α and β are Cartesian polarization, e β (s|υ ) indicates the normalized vibrational eigenvector, and Z * αβ indicates the effective charge tensor. To assess how the frequencies obtained depend on the computational setting, the performance of four different density functionals (PBE, Local Density Approximation (LDA), Perdew-Wang (PW91), Perdew-Burke-Ernzerhof revised for solids (PBESOL)), cut-off (300, 400, 500, 600, and 700 eV), choice of U (3, 4, 5, 6, and 7 eV), and the inclusion of dipole corrections were tested (see Supplementary Tables S1-S4 and Figures S7-S10). Although the numerical values are affected by computational settings, the trends between the different orientations are maintained. Dispersion effects were evaluated for the heterolytic path for the intermediates and TS1 (the latter as a single-point calculation) by means of dispersion corrections (Grimme D3) (zero) [80] and the results are displayed in Supplementary Table S7 and Figure S12. Dispersion corrections were found to slightly stabilize adsorbates with respect to the non-corrected calculation and significantly decrease the barriers. However, they did not significantly alter the trends of the terminations, nor the vibrational frequencies (Supplementary Table S7 and Figure S12).
No dipole correction was used to account for the asymmetry of the slabs in the perpendicular direction. As our work is mainly based on the comparison of terminations, and as all of them should be affected in a similar manner by the spurious dipole, we do not expect it to have a significant impact on the conclusions. As can be seen in Supplementary Table S6 and Figure S11, the inclusion of dipole corrections did not have a significant effect in the vibrational frequencies.

Slab Model
We optimized the bulk TiO 2 rutile unit cell obtaining values of a = b = 4.661 Å and c = 2.961 Å, in agreement with experimental parameters of a = 4.593 Å and c = 2.958 Å [23]. The calculated lattice parameters for bulk rutile TiO 2 were overestimated by 1.46% for a and only 0.10% for c with respect to the experimental value, and the optimized values were used to build the slab models.
The four rutile TiO 2 surface (001), (100), (110), and (101) stoichiometric terminations are represented in Figure 1 and the main structural parameters are reported in Table 1. As we can see, the facets (001) and (110) are roughly flat, while (100) and (101) facets are uneven. On the surfaces, the coordination number of titanium sites vary from 4 to 6(001) has only Ti 4C ; (101) possesses Ti 4C and Ti 5C ; (100) has only Ti 5C ; and (110) has Ti 5C and Ti 6C . Regarding oxygen, the surface coordination varies from two-to three-fold-(001) has only O 2C , while the other three exhibit O 2C and O 3C . The surface energy E surf , calculated as the difference in energy between the slab and the bulk divided by twice the area, follows the trend of coordination-the lower the surface atomic coordination, the higher the surface energy. Thus, (001), where Ti and O are poorly coordinated, shows E surf 1.30 J nm −2 , whereas (110), where the atoms are more coordinated, shows E surf 0.55 J nm −2 .

H 2 Dissociation
Firstly we investigated the heterolytic pathway for H 2 dissociation on the four selected rutile TiO 2 selected. In the first step, the H 2 molecule physisorbs on the surface forming the adduct H 2 *. Then, the heterolytic H 2 dissociation takes place between the Ti site and a neighboring O atom through a transition structure (TS1), generating a pair of O−H and Ti−H bonds (H + -H − species). The second step involves the transfer of the hydride (H − ) on the Ti site to a nearby O, leading to 2 O-H hydroxyl groups (H + -H + ) and a two-electron transfer to surface titanium sites that become reduced. The transition state associated with this step is labeled as TS2. The reaction pathway involving these two steps is schematized in Figure 2, and the calculated energies are reported in Table 2. The energy profile of H 2 dissociation over the four surfaces is shown in Figure 3. found to slightly stabilize adsorbates with respect to the non-corrected calculation and significantly decrease the barriers. However, they did not significantly alter the trends of the terminations, nor the vibrational frequencies (Supplementary Table S7 and Figure S12). No dipole correction was used to account for the asymmetry of the slabs in the perpendicular direction. As our work is mainly based on the comparison of terminations, and as all of them should be affected in a similar manner by the spurious dipole, we do not expect it to have a significant impact on the conclusions. As can be seen in Supplementary Table S6 and Figure S11, the inclusion of dipole corrections did not have a significant effect in the vibrational frequencies.

Slab Model
We optimized the bulk TiO2 rutile unit cell obtaining values of a = b = 4.661 Å and c = 2.961 Å, in agreement with experimental parameters of a = 4.593 Å and c = 2.958 Å [23]. The calculated lattice parameters for bulk rutile TiO2 were overestimated by 1.46% for a and only 0.10% for c with respect to the experimental value, and the optimized values were used to build the slab models.
The four rutile TiO2 surface (001), (100), (110), and (101) stoichiometric terminations are represented in Figure 1 and the main structural parameters are reported in Table 1. As we can see, the facets (001) and (110) are roughly flat, while (100) and (101) facets are uneven. On the surfaces, the coordination number of titanium sites vary from 4 to 6(001) has only Ti4C; (101) possesses Ti4C and Ti5C; (100) has only Ti5C; and (110) has Ti5C and Ti6C. Regarding oxygen, the surface coordination varies from two-to three-fold-(001) has only O2C, while the other three exhibit O2C and O3C. The surface energy Esurf, calculated as the difference in energy between the slab and the bulk divided by twice the area, follows the trend of coordination-the lower the surface atomic coordination, the higher the surface energy. Thus, (001), where Ti and O are poorly coordinated, shows Esurf 1.30 J nm −2 , whereas (110), where the atoms are more coordinated, shows Esurf 0.55 J nm −2 .

H2 Dissociation
Firstly we investigated the heterolytic pathway for H2 dissociation on the four selected rutile TiO2 selected. In the first step, the H2 molecule physisorbs on the surface forming the adduct H2*. Then, the heterolytic H2 dissociation takes place between the Ti site and a neighboring O atom through a transition structure (TS1), generating a pair of O−H and Ti−H bonds (H + -H − species). The second step involves the transfer of the hydride (H − ) on the Ti site to a nearby O, leading to 2 O-H hydroxyl groups (H + -H + ) and a two-electron transfer to surface titanium sites that become reduced. The transition state associated with this step is labeled as TS2. The reaction pathway involving these two steps is schematized in Figure 2, and the calculated energies are reported in Table 2. The energy profile of H2 dissociation over the four surfaces is shown in Figure 3.     Here, all adsorption energies are referred to the energy of physisorbed TiO2-H2. The path for the TiO2 (001) surface is illustrated, and those corresponding to the other three terminations are provided in Supplementary Figures S1-S3. Several pathways were considered involving different surface sites. For the heterolytic step, on (001) there is a unique possible pathway with only one kind of O2c site and one Ti4c site on the surface. On (100), besides the pathway reported in Supplementary Figure S1, there is also one additional combination of Ti5C and O2C sites (Supplementary Figure S6a), in which the direction of the OH bond is almost perpendicular to the direction of the TiH bond, which makes this combination less stable than the one selected. For (101), two other possible structures involving Ti5C and O2C (Supplementary Figure S6b) and Ti4C and O3C (Supplementary Figure S6c) resulted in less stable systems than the one retained. The model structures retained are stabilized as a consequence of the saturation of poorly coordinated sites of Ti4C, Ti5C, and O2C upon hydrogenation, and in some cases the formation of hydrogen bonds.  (001) and (101) terminations (−0.08 eV), whereas it is slightly endothermic for the (110) by 0.12 eV, and for the (100) by 0.68 eV. Hydrogen bonds between TiH and OH species form in all the terminations except (101). Whereas the terminations showing the poorest coordination exhibit the most exothermic adsorption energy for the (H + , H − ) intermediate, the most highly coordinated slabs show less exothermic values. However, the most highly coordinated (110) slab exhibits a significantly lower adsorption energy than (100). The activation barriers of heterolytic H 2 dissociation on the four TiO 2 surfaces follow the trend (110) 0.50 eV < (001) 0.56 eV < (101) 0.79 eV < (100) 1.08 eV. As for the adsorption energy, a trend appears between coordination and kinetic barriers for (001), (101), and (100), whereas (110) presents lower values than expected (its higher coordination should lead to the most endothermic values). Our results are consistent with previous studies. For the (001) surface, our activation energy (0.56 eV) is consistent with the one reported previously (0.68 eV) [14]. The difference comes from the use of a different Hubbard parameter (U = 7 eV) and unit cell (2 × 1). Our activation energy for the (110) surface, 0.50 eV, is larger than the 0.37 eV reported for a much narrower slab (3-TiO 2 -layer thick slab model [34]), highlighting the important role of slab thickness in the construction of a model.
According to our results, the (H + , H − ) intermediate is more likely to be formed on (001) and (101) terminations, however the poor stability and the low barriers could induce the inverse reaction, i.e., the recombination and desorption as H 2 (see below). These results suggest that the rutile TiO 2 (001) exhibits the most likely H 2 heterolytic dissociation path in the series, with the lowest activation energy and a slight stabilization of the product. This specific reactivity could be associated with the low coordination of the surface titanium site-the four-fold coordinated Ti site in the (001) termination stabilizes the hydride Ti-H species to increase the number of neighbors. In the transition structure 1 (TS1) displayed in Figure 3 However, the most highly coordinated (110) slab exhibits a significantly lower adsorption energy than (100). The activation barriers of heterolytic H2 dissociation on the four TiO2 surfaces follow the trend (110) 0.50 eV < (001) 0.56 eV < (101) 0.79 eV < (100) 1.08 eV. As for the adsorption energy, a trend appears between coordination and kinetic barriers for (001), (101), and (100), whereas (110) presents lower values than expected (its higher coordination should lead to the most endothermic values). Our results are consistent with previous studies. For the (001) surface, our activation energy (0.56 eV) is consistent with the one reported previously (0.68 eV) [14]. The difference comes from the use of a different Hubbard parameter (U = 7 eV) and unit cell (2 × 1). Our activation energy for the (110) surface, 0.50 eV, is larger than the 0.37 eV reported for a much narrower slab (3-TiO2-layer thick slab model [34]), highlighting the important role of slab thickness in the construction of a model.
According to our results, the (H + , H − ) intermediate is more likely to be formed on (001) and (101) terminations, however the poor stability and the low barriers could induce the inverse reaction, i.e., the recombination and desorption as H2 (see below). These results suggest that the rutile TiO2 (001) exhibits the most likely H2 heterolytic dissociation path in the series, with the lowest activation energy and a slight stabilization of the product. This specific reactivity could be associated with the low coordination of the surface titanium site-the four-fold coordinated Ti site in the (001) termination stabilizes the hydride Ti-H species to increase the number of neighbors. In the transition structure 1 (TS1) displayed in Figure 3 (Figure 4). This is consistent with a moderate polarization of the H2 moiety, as shown in the Bader analysis discussed below. In this TS1 structure, the four atoms involved Ti-H…H-O are coplanar.   (100), the final products involve a rearrangement of the surface bonds-a Ti-O-Ti breaks to form a Ti-OH moiety. The activation barriers are significantly higher than for the first step, ranging from 1.30 to 1.80 eV. These results indicate an unfavorable evolution to the homolytic product from the heterolytic intermediate. Thus, the hydride TiH species could be kinetically stabilized on TiO 2 surfaces with a possible recombination to regenerate and desorb H 2 at low temperatures, whereas the reduction step would require much higher energies to occur. Nevertheless, the most thermodynamically stable product is found for step 2 and involves the presence of two hydroxyl groups and two Ti 3+ sites; the latter originate from the electron transfer from the hydride to two titanium sites. This transfer results in open-shell systems that can be characterized by the presence of two unpaired electrons.
We have looked for correlations between adsorption energy, barrier heights, and geometry (TiH, HH, and OH distances), as well as Bader charges in TSs, and our results indicate no clear relationship. This is very interesting, as for CeO 2 those correlations do appear [81]. This might point to an ionocovalent character of Ti-O bond compared to the more ionic Ce-O bond, which would facilitate the formation of the H + -H − ionic pair, or to the important role of the local topology in stabilizing intermediates and transition structures. As a general trend, the activation barriers seem related to the coordination numbers of Ti and O on the surfaces, with the (110) termination behaving in a different way than is expected from its highly coordinated surface sites.

Electronic Structure
In order to characterize in more detail the electronic structure of the structures involved in the hydrogenation mechanisms, we have computed the density of states (DOS); Figure 5 Figure S4. At the bottom of the plot a hydrogen molecular band appears as a sharp narrow peak in the valence region due to H 2 physical adsorption. In TS1, we observe a splitting in two bands associated with H + and H − species that overlap with the slab levels. For the product of heterolytic dissociation (H + -H − ), the H − band is the highest occupied energy level, with a sharp peak at the Fermi level. For the TS2 of subsequent H transfer from Ti to nearby O, there still exists one H + band and one H − band, but the intensity decreases. The existence of wide, weak peaks in the gap indicates an early reduction of the Ti site in TS2 on (100) and (110) facets, while no corresponding peak appears on (001) and (101). For the H transfer process product (H + -H + ) species, we observe the H + levels corresponding to OH groups in O-H bonds in the valence band, which appear as two distinct peaks if they correspond to inequivalent hydroxyl groups. Also, Ti states appear in the gap below the Fermi level, indicating the reduction of the Ti sites. This is consistent with the picture of the spin density plots (Figure 6 and Supplementary Figure S5), indicating that the unpaired electrons from the hydride transfer are trapped by two Ti ions that get reduced, confirming the nature of Ti 3+ sites. Note that the approach used in the present work does not allow one to state unambiguously which Ti sites are reduced-it only confirms qualitatively that two distinct Ti sites are involved.
Bader charge analysis [82] was carried out to complement the characterization of the electronic structure of the systems studied. In Table 3, we can follow the electronic charges during the two processes, whereas Table 4 shows the Bader analysis for the spin density. In step 1, the adsorbed hydrogen species shows a slightly polarized H-H bond. In TS1, the H-H bond is more polarized,  Table 4) two Ti are reduced for every facet. For the TS2, one of the Ti on the surface partially decreases its positive charge, indicating partial reduction. Finally, in the H + -H + species the two H are characterized as protons, whereas two Ti sites decrease their positive charge, indicating that they host the reduction electrons, and the integrated spin density varies from 0.90 to 1.05 |e| (See Table 4). It is worth stating that the O site involved in the H transfer process also contains a small amount of unpaired electrons of 0.24 |e| for TS2 of both (001) and (101) facets, as can be seen in Figure 6 for the (001) case.    Bader charge analysis [82] was carried out to complement the characterization of the electronic structure of the systems studied. In Table 3, we can follow the electronic charges during the two processes, whereas Table 4 shows the Bader analysis for the spin density. In step 1, the adsorbed hydrogen species shows a slightly polarized H-H bond. In TS1, the H-H bond is more polarized,

Effect of the Hubbard Correction U
The values without U correction were considered to analyze the effect of U in the energetic profile, which is reported in Table 2. The profile is similar to the one obtained for U = 4 eV (see Figure 7 for the (001) case and Supplementary Materials for the others). In step 1, the heterolytic dissociation leads to (H + -H − ) products stable at 0.15 eV (001), 0.28 eV (101), 0.98 eV (101), and 0.50 eV (110), and barriers of 0.63 eV, 1.10 eV, and 1.15 eV, 0.70 eV, respectively, which is~0.20 eV higher in energy than for the U = 4 eV case (Table 2, Figure 3). The increase in the values is significantly higher in step 2, where the (H + -H + ) product is higher in energy by~0.60 eV in the absence of U correction, and is associated with the stabilization of the localized solution favored by the U = 4 eV term with respect to the U = 0 eV case. In general, the activation barriers are not significantly affected by the U value, with the exception of (110) and (100)

Vibrational Spectrum
We computed the vibration frequency and IR spectra of H2 heterolytic dissociation products (H + -H − ) for the four terminations. No scaling factor was applied. The vibration modes are shown in Table 5 and Figure 8, and present three main regions: Ti-H and O-H bending modes at low frequencies (below 1000 cm −1 ); the vibration frequencies of Ti-H lie in the range of 1500-1800 cm −1 ; the stretching OH modes are characterized by higher frequencies (between 2900 cm −1 and 3800 cm −1 ). Ti-H stretching modes of the four species are seen in the calculated spectrum at 1644 cm −1 (001), 1768 cm −1 (100), 1653 cm −1 (110), and 1577 cm −1 (101), corresponding to the expected Ti-H IR spectral region (around 1600 cm −1 ) [83].The hydrides of (001) and (101) facets are Ti4C-H, while they are Ti5c-H on (100) and (110) surfaces, as displayed in Figure 9. Previous studies using electron-stimulated desorption (ESD) [84] and low-energy ion scattering (LEIS) [85] reported that the annealed TiO2 surface is compensated by H, which is bonded in the Ti-H as well as O-H with bridging O or a subsurface, but no specific frequencies were provided. Recently, Yan et.al indicated the formation of Ti-H species on the P25 TiO2 surface [86]. Table 5. Computed IR wavenumbers (cm −1 ) and intensities (in brackets) of Ti−H and O−H stretching modes of (H + , H − ) species for the four terminations studied.

Vibrational Spectrum
We computed the vibration frequency and IR spectra of H 2 heterolytic dissociation products (H + -H − ) for the four terminations. No scaling factor was applied. The vibration modes are shown in Table 5 and Figure 8, and present three main regions: Ti−H and O−H bending modes at low frequencies (below 1000 cm −1 ); the vibration frequencies of Ti−H lie in the range of 1500-1800 cm −1 ; the stretching OH modes are characterized by higher frequencies (between 2900 cm −1 and 3800 cm −1 ). Ti−H stretching modes of the four species are seen in the calculated spectrum at 1644 cm −1 (001), 1768 cm −1 (100), 1653 cm −1 (110), and 1577 cm −1 (101), corresponding to the expected Ti−H IR spectral region (around 1600 cm −1 ) [83].The hydrides of (001) and (101) facets are Ti 4C -H, while they are Ti 5c -H on (100) and (110) surfaces, as displayed in Figure 9. Previous studies using electron-stimulated desorption (ESD) [84] and low-energy ion scattering (LEIS) [85] reported that the annealed TiO 2 surface is compensated by H, which is bonded in the Ti−H as well as O−H with bridging O or a subsurface, but no specific frequencies were provided. Recently, Yan et.al indicated the formation of Ti-H species on the P25 TiO 2 surface [86].

Vibrational Spectrum
We computed the vibration frequency and IR spectra of H2 heterolytic dissociation products (H + -H − ) for the four terminations. No scaling factor was applied. The vibration modes are shown in Table 5 and Figure 8, and present three main regions: Ti-H and O-H bending modes at low frequencies (below 1000 cm −1 ); the vibration frequencies of Ti-H lie in the range of 1500-1800 cm −1 ; the stretching OH modes are characterized by higher frequencies (between 2900 cm −1 and 3800 cm −1 ). Ti-H stretching modes of the four species are seen in the calculated spectrum at 1644 cm −1 (001), 1768 cm −1 (100), 1653 cm −1 (110), and 1577 cm −1 (101), corresponding to the expected Ti-H IR spectral region (around 1600 cm −1 ) [83].The hydrides of (001) and (101) facets are Ti4C-H, while they are Ti5c-H on (100) and (110) surfaces, as displayed in Figure 9. Previous studies using electron-stimulated desorption (ESD) [84] and low-energy ion scattering (LEIS) [85] reported that the annealed TiO2 surface is compensated by H, which is bonded in the Ti-H as well as O-H with bridging O or a subsurface, but no specific frequencies were provided. Recently, Yan et.al indicated the formation of Ti-H species on the P25 TiO2 surface [86].    Table S5), showing that the vibrations can be greatly influenced by the nature of the site, the surface topology, the presence of defects and coverage [87], as well as the polymorph [88]. In the present work, we perform an analysis of OH vibrations for the four terminations considered (see Table 5, Figure 9) as a guide for qualitative assignment. It is found that the OH stretching vibrations are different between these four facets.  [89]. For the TiO2 (110) surface (3606.59 cm −1 ) the vibration corresponds to Ti6C-O2CH, which is lower than in a previous theoretical study by Wöll (3700 cm −1 ) [90]. Note that the model used in the work of Wöll et al. involves a hydroxyl perpendicular to the slab, whereas in our work the hydroxyl is tilted. Other experimental works report 3665 and 3690 cm −1 measured by High-Resolution Electron Energy Loss Spectroscopy (HREELS) [91,92], and 3711 cm −1 by IR [93] on systems obtained by H2O adsorption on a clean single-crystal TiO2 surface. As a general remark, the lack of experimental data in well-controlled structures and conditions make an assessment of the vibrational spectra of surface hydroxyl and hydride species difficult, although several trends can be observed. First, the vibrations are dependent on the surface topology due to specific local chemical environments. Second, the coordination of oxygen and titanium sites seems to play a role, as well as hydrogen bonds formed between TiH/OH pairs and neighboring O sites. Overall, our results are consistent with previous experimental and theoretical data published in the literature and provide a set of spectra to stimulate the search of TiH/OH species on different rutile terminations.

The H2 Recombination-Desorption Reaction
We studied the energy barriers for hydride TiH/OH species recombination to regenerate and desorb H2 on four facets ( Figure 10, Table 2). The corresponding barrier for that process,   Table S5), showing that the vibrations can be greatly influenced by the nature of the site, the surface topology, the presence of defects and coverage [87], as well as the polymorph [88]. In the present work, we perform an analysis of OH vibrations for the four terminations considered (see Table 5, Figure 9) as a guide for qualitative assignment. It is found that the OH stretching vibrations are different between these four facets. The calculated IR results of TiO 2 (100) surface (2976.54 cm −1 ) correspond to a Ti 5C -O 3C H exhibiting an H bond with one O site nearby. An experimental value of 3550 cm −1 for OH on (100) [89] was reported for the adsorption of water on the surface, most likely assigned to terminal hydroxyl groups. Our value is consistent with a higher coordination of the hydroxyl group (three-fold in our case), as well as with the presence of a hydrogen bond, both blue-shifting the vibration with respect to the experimental value. For the TiO 2 (101) surface (3622.37 cm −1 ) it corresponds to a Ti 4C -O 2C H. Experimentally, the OH stretching vibrations from water adsorption are observed at 3680 and 3610 cm −1 [89]. For the TiO 2 (110) surface (3606.59 cm −1 ) the vibration corresponds to Ti 6C -O 2C H, which is lower than in a previous theoretical study by Wöll (3700 cm −1 ) [90]. Note that the model used in the work of Wöll et al. involves a hydroxyl perpendicular to the slab, whereas in our work the hydroxyl is tilted. Other experimental works report 3665 and 3690 cm −1 measured by High-Resolution Electron Energy Loss Spectroscopy (HREELS) [91,92], and 3711 cm −1 by IR [93] on systems obtained by H 2 O adsorption on a clean single-crystal TiO 2 surface.
As a general remark, the lack of experimental data in well-controlled structures and conditions make an assessment of the vibrational spectra of surface hydroxyl and hydride species difficult, although several trends can be observed. First, the vibrations are dependent on the surface topology due to specific local chemical environments. Second, the coordination of oxygen and titanium sites seems to play a role, as well as hydrogen bonds formed between TiH/OH pairs and neighboring O sites. Overall, our results are consistent with previous experimental and theoretical data published in the literature and provide a set of spectra to stimulate the search of TiH/OH species on different rutile terminations.

The H 2 Recombination-Desorption Reaction
We studied the energy barriers for hydride TiH/OH species recombination to regenerate and desorb H 2 on four facets ( Figure 10, Table 2). The corresponding barrier for that process, E back act , requires 0.64 eV for (001), 0.87 eV for (101), 0.38 eV for the (110), and 0.40 eV for the (100) slabs. The backward activation energies for the facets (001) and (101) are larger than those found for facets (110) and (100), probably due to the higher stability of the (H + -H − ) species. Contrary to the dissociation process, the desorption of H 2 is slightly endothermic for the (001) and (101) terminations (0.08 eV), whereas it is exothermic for the (110) by −0.12 eV, and for the (100) by −0.

Zero-Point Energy Correction and Effect of Temperature
The energy profiles with Zero Point Energy (ZPE) correction are also studied together with the Gibbs free energies for T = 298 K ( Figure 11). With ZPE correction, the energy for these two steps increases, while it does not affect the kinetic barriers. Temperature has almost no effect on this reaction profile. As a final remark, many other factors may have a deep influence on the behavior of TiO2 regarding hydrogenation-the presence of surface and subsurface defects [94], the nature of the bulk phase [95], nanostructuring [96,97], interfacial water [98], or reduction [99]. More fundamental works to elucidate the structure of hydrogenated surfaces are needed to build a robust scenario for the complex behavior observed [100].

Conclusions
The mechanisms of H2 dissociation on four different rutile TiO2 facets by means of density functional theory (PBE+U) calculations have been investigated. The results showed that the topology of the surface has a moderate effect on H2 dissociation on TiO2 kinetically and also

Zero-Point Energy Correction and Effect of Temperature
The energy profiles with Zero Point Energy (ZPE) correction are also studied together with the Gibbs free energies for T = 298 K ( Figure 11). With ZPE correction, the energy for these two steps increases, while it does not affect the kinetic barriers. Temperature has almost no effect on this reaction profile.

Zero-Point Energy Correction and Effect of Temperature
The energy profiles with Zero Point Energy (ZPE) correction are also studied together with the Gibbs free energies for T = 298 K ( Figure 11). With ZPE correction, the energy for these two steps increases, while it does not affect the kinetic barriers. Temperature has almost no effect on this reaction profile. As a final remark, many other factors may have a deep influence on the behavior of TiO2 regarding hydrogenation-the presence of surface and subsurface defects [94], the nature of the bulk phase [95], nanostructuring [96,97], interfacial water [98], or reduction [99]. More fundamental works to elucidate the structure of hydrogenated surfaces are needed to build a robust scenario for the complex behavior observed [100].

Conclusions
The mechanisms of H2 dissociation on four different rutile TiO2 facets by means of density functional theory (PBE+U) calculations have been investigated. The results showed that the topology of the surface has a moderate effect on H2 dissociation on TiO2 kinetically and also As a final remark, many other factors may have a deep influence on the behavior of TiO 2 regarding hydrogenation-the presence of surface and subsurface defects [94], the nature of the bulk phase [95], nanostructuring [96,97], interfacial water [98], or reduction [99]. More fundamental works to elucidate the structure of hydrogenated surfaces are needed to build a robust scenario for the complex behavior observed [100].

Conclusions
The mechanisms of H 2 dissociation on four different rutile TiO 2 facets by means of density functional theory (PBE+U) calculations have been investigated. The results showed that the topology of the surface has a moderate effect on H 2 dissociation on TiO 2 kinetically and also thermodynamically. We found that for all four surfaces, the heterolytic dissociation pathway towards hydride-hydroxyl surface pairs is kinetically more favorable than the H transfer process towards substrate reduction, although the reduction product, with only surface hydroxyl groups, is thermodynamically more favorable. On (110) and (100), the hydride-hydroxyl pair formed can recombine and desorb as molecular dihydrogen, whereas the (001), and to a lesser extent (101), stabilize the hydride-hydroxyl pair. The energetics of the reaction seems related to the coordination numbers of Ti and O on the surfaces, although (110) shows a specific behavior. No clear trend relating adsorption energies and barriers with local geometry or charges was found. The electronic structure analysis allows characterization of charge and electron transfers. The IR spectra of the (H + -H − ) pair species were also computed indicating the vibrational region of Ti-H species on TiO 2 facets in the range of 1550-1750 cm −1 . The frequencies are found to depend on the facet exposed and could be used as a qualitative guideline to identify them experimentally.