Breaking H2 with CeO2: Effect of Surface Termination.

The ability of ceria to break H2 in the absence of noble metals has prompted a number of studies because of its potential applications in many technological fields. Most of the theoretical works reported in the literature are focused on the most stable (111) termination. However, recently, the possibility of stabilizing ceria particles with selected terminations has opened new avenues to explore. In the present paper, we investigate the role of termination in H2 dissociation on stoichiometric ceria. We model (111)-, (110)-, and (100)-terminated slabs together with the stepped (221) and (331) surfaces. Our results support a dissociation mechanism proceeding via the formation of a hydride/hydroxyl CeH/OH intermediate. Both the stability of such an intermediate and the activation energy depend critically on the termination, the (100)-terminated surfaces being the most reactive: the activation energy is 0.16 eV, and the CeH/OH intermediate is stable by -0.64 eV for the (100) slab, whereas the (111) slab presents 0.75 and 0.74 eV, respectively. We provide structural, energetic, electronic, and spectroscopic data, as well as chemical descriptors correlating structure, energy, and reactivity, to guide in the theoretical and experimental characterization of the Ce-H surface intermediate.

Introduction

Ceria has attracted much attention in the last years because of its numerous technological application fields such as heterogeneous catalysis,[1,2] where it can be used as a catalyst itself[3−6] or as a support,[7−13] treatment of toxic gases and pollutants,[14,15] solid oxide fuel cells,[16,17] oxygen sensors,[18−20] and biomedicine.[21,22] Pure ceria has been successfully used in alkyne semihydrogenation reactions[5,23−25] with high activity and selectivity to the alkene products. The unexpected ability of ceria to dissociate hydrogen opens new directions for the use of this promising material, in particular in the field of heterogeneous catalysis where the absence of noble metal particles involves tremendous economic advantages. These last years, H2 dissociation on ceria has been of great interest both experimentally[24−31] and theoretically.[24,25,27,32−37] If it is well-established that the product of hydrogenation in ceria leads to stable hydroxylated surfaces,[30,34,36,38] then the mechanism was found to proceed through a hydride intermediate: H2 dissociates on ceria according to a heterolytic pathway (forming a hydride/hydroxyl pair), followed by the transfer of a H atom that finally yields the homolytic product (forming two hydroxyl groups and reducing two cerium atoms). In particular, García-Melchor and López[34] and Fernández-Torre et al.[36] find that the H–H bond break on (111)-CeO2 proceeds via the formation of a metastable intermediate Ce–H/Ce–OH, which is ∼0.7 eV higher in energy compared to H2 molecular adsorption. Therefore, understanding the dissociation mechanism of H2 on ceria becomes crucial not only to rationalize the elementary steps involved, but also for the development of materials with selected properties such as heterogeneous catalysts. Heterolytic bond break pathways are commonly associated in the literature to irreducible oxides,[39,40] and the fact that they drive the H2 dissociation in ceria opens new avenues to explore.

Recently, some of us[37] have conducted a detailed analysis of the dissociation mechanism on ceria and gallia-promoted ceria catalysts by a multitechnique approach involving X-ray photoelectron spectroscopy (XPS), infrared (IR), nuclear magnetic resonance (NMR), microkinetic modeling and density functional theory (DFT). Our results show that the formation of an intermediate CeH/OH pair is consistent with the activation energy deduced from the experimental measurements, whereas a direct dissociation to the final hydroxylated product is not supported by the data. In the mechanism proposed, the formation of Ce–H hydride species is the key step. Such species were not observed in ceria by IR spectroscopy;[25,28,37] on the contrary, Ga–H species are clearly identified in the doped materials in the same set of experiments[37] and also in other experiments[41,42] and are associated to lower activation barriers for the hydrogenation reaction. H2 dissociation was also studied on ceria-supported gold nanoparticles, where it has been shown that H2 dissociates according to a heterolytic pathway, leading to the formation of Au–H and O–H species.[43]

In all the theoretical works cited above, the most stable (111) termination of ceria is considered as model for the mechanistic studies of H2 dissociation. In standard conditions and under thermodynamic equilibrium, ceria particles will expose mostly the facets corresponding to the most stable surfaces, whereas the other terminations will count for a small fraction of the surface particle. Thus, octahedral nanoparticles expose mostly (111) facets that can be truncated depending on the external conditions.[44−46] However, other terminations of ceria are experimentally accessible, such as (100) which is predominant in nanocubes,[29,47−55] (110) in nanorods,[52−58] and the most common (111) in nano-octahedra[29,52−55] ceria particles. In the present study, the roles of surface termination as well as surface topology in the H2 bond break mechanism were investigated, showing in particular that the Ce–H species can be stabilized, leading to a decrease in the activation barrier of H2 break. In particular, we have shown that a careful selection of the surface termination can decrease up to 5 times the activation barrier. This study focuses on stoichiometric surfaces, and the reduced surfaces will be addressed in future works. The (100), (110), and (111) surfaces were investigated, together with the stepped (221) and (331) ones, which were also reported to be stable.[59] Our results suggest that the (100)-terminated surfaces are the most efficient to break the H–H bond, with a barrier as low as 0.16 eV. Besides, the heterolytic dissociation product (O–H and Ce–H species) on these surfaces is found to be thermodynamically stable (between −0.20 and −0.64 eV), and we provide computed IR spectra as well as temperature effects (Gibbs free energies) that may guide the experimentalists in the search of the characterization of this important intermediate.

The present paper is organized as follows. The results for the bare slabs and their interaction with H2 are presented, followed by the electronic structure analysis and the effect of temperature. The use of reactivity descriptors is discussed. The conclusions raised and the methodology employed are presented at the end of the paper.

Results and Discussion

Bare Slabs

The seven surface models used are depicted in Figure . From a structural point of view, (100)-Ce, (110), and (111) are terminated by threefold-coordinated oxygens, and the (100)-O-(a)/(b) slabs show twofold-coordinated oxygens. The coordination number of surface cerium sites in (100)-Ce is fourfold; in (100)-O-(a)/(b) and (110), it is sixfold, and in (111), it is sevenfold-coordinated. The step sites are characterized by threefold-coordinated oxygen and sixfold-coordinated cerium (Figure and Table ).

Figure 1

Side view of the structure for the different terminations studied. Cerium and oxygen atoms are depicted by yellow and red spheres, respectively. For (221) and (331) surfaces, the step and (111) “type” sites are shown by a blue square and a purple rectangle, respectively.

Table 1

Size, Composition (CeO2 Units), Number of Atomic Layers and Area (nm2) of the Supercell, Coordination Number, and Number of Surface Sites per Supercell of Surface Sites as well as Surface Energy (J/m2) for the Slab Models

	(100)-Ce	(100)-O-(a)/(b)	(110)	(111)	(221)	(331)
supercell	2 × 2	2 × 2	2 × 2	3 × 3	1 × 2	1 × 1
composition (CeO2 units)	72	72	64	36	52	40
atomic layers (frozen/relaxed)	19 (10/9)	19 (10/9)	8 (4/4)	12 (6/6)	27 (14/13)	30 (15/15)
area (nm2)	1.21	1.21	1.71	1.18	1.82	1.32
coordination number of surface sites	O(3)	O(2)	O(3)	O(3)	O(3)	O(3)
	Ce(4)	Ce(6)	Ce(6)	Ce(7)	Ce(6; 7)	Ce(6; 7)
surface site number	O(16)	O(8)	O(12)	O(9)	O(16)	O(12)
	Ce(4)	Ce(8)	Ce(8)	Ce(9)	Ce(4; 8)	Ce(4; 4)
γhkl (J/nm2)	1.83	1.48/1.64	1.09	0.73	0.88	0.95

The surface energies were evaluated according towhere Eslab, is the energy of the (hkl) supercell, N is the number of CeO2 bulk unit in the supercell, Ebulk is the reference energy of one bulk unit, and A is the surface area of the supercell.

H2 Dissociation

A mechanism involving different steps is schematized in Figure , and the calculated energies are reported in Table . First, the H2 molecule physisorbs at an oxygen site (R1) and is preactivated, with a small polarization, as can be seen in Table . In a second step, the H–H bond is broken according to a heterolytic path leading to a metastable state, labeled MS, and achieved after a process (R2′), that evolves spontaneously, that is, without an energetic barrier, to the product of the heterolytic dissociation through the process (R2″). The product is a pair formed by a hydride bound to the Ce site and a proton bound to O, labeled (H+, H–). For the (100) surfaces, H2 dissociates directly without MS, via the process (R2). Finally, because of the strong exothermicity of the formation of hydroxyls, we consider the product of the homolytic dissociation with the diffusion of one H atom to an oxygen site to form two hydroxyl groups—process (R3).

Figure 2

Schematic representations of the different steps involved in the heterolytic dissociation pathway. The species are labeled as follows: H2 physisorbed as H2*, metastable state as MS, hydride–proton pair from the heterolytic product as (H+, H–), and the homolytic product as (OH, OH).

Table 2

Adsorption Energies, Reaction Energies (Erea), Forward Activation Energies (Eactfor), and Backward Activation Energies (Eactback) in Reference to H2 Dissociation Pathwaya

	(100)-Ce	(100)-O-(a)	(100)-O-(b)	(110)	(111)	(221)	(331)
H2*	–0.15 (0.35)	–0.09 (0.41)	–0.14 (0.36)	–0.08 (0.42)	–0.08 (0.40)	0.08 (0.40)	–0.08 (0.40)
TS	0.20 (0.64)	0.30 (0.75)	0.02 (0.46)	0.45 (0.89)	0.67 (1.10)	0.50 (0.93)	0.53 (0.96)
MS				0.43 (0.92)	0.66 (1.11)	0.47 (0.95)	0.51 (0.99)
(H+, H–)	–0.64 (−0.13)	–0.27 (0.27)	–0.20 (0.30)	0.34 (0.86)		0.46 (0.96)	0.48 (0.98)
(OH, OH)	–2.35 (−1.67)	–3.47 (/)	–3.85 (/)	–3.19 (−2.51)	–2.18 (−1.51)	–3.06 (−2.38)	–2.86 (−2.17)
Erea	–0.49 (−0.48)	–0.18 (−0.14)	–0.06 (−0.06)	0.42 (0.44)	0.74 (0.71)	0.54 (0.56)	0.56 (0.58)
Eactforw	0.35 (0.29)	0.39 (0.34)	0.16 (0.10)	0.53 (0.47)	0.75 (0.70)	0.58 (0.53)	0.61 (0.56)
Eactback	0.84 (0.77)	0.57 (0.48)	0.22 (0.16)	0.11 (0.03)	0.01 (−0.01)	0.04 (−0.03)	0.05 (−0.02)

All adsorption energies are referenced to the total energy of H2(g) and the stoichiometric CeO2 surface. Asterisk indicates the adsorbed molecule (not in gas phase). Values in bracket are Gibbs free energies for T = 300 K and PH = 1 bar.

Table 4

Bader Charges (in |e|) of Selected Atoms Involved in the Dissociation Process for Slab, H2*, TS, and (H+, H–) Species (For the (111) Surface, Bader Charges were Given for MS Instead of (H+, H−))

	slab	H2*	TS	(H+, H–)
	qCe/qO	qHδ+/qHδ−/qCe/qO
(100)-Ce	+2.05/–1.09	+0.13/–0.15/+2.07/–1.11	+0.39/–0.49/+2.14/–1.13	+1.00/–0.68/+2.24/–1.65
(100)-O-(a)	+2.34/–1.12	+0.02/–0.02/+2.34/–1.12	+0.36/–0.39/+2.34/–1.14	+1.00/–0.46/+2.31/–1.67
(100)-O-(b)	+2.32/–1.11	+0.07/–0.08/+2.32/–1.11	+0.34/–0.41/+2.33/–1.13	+1.00/–0.56/+2.32/–1.67
(110)	+2.32/–1.16	+0.02/–0.01/+2.32/–1.15	+0.43/–0.46/+2.32/–1.20	+1.00/–0.51/+2.31/–1.68
(111)	+2.37/–1.18	+0.03/–0.03/+2.40/–1.18	+1.00/–0.44/+2.35/–1.76	+1.00/–0.43/+2.33/–1.77
(221)	+2.33/–1.17	+0.02/–0.01/+2.33/–1.17	+0.45/–0.44/+2.31/–1.22	+1.00/–0.48/+2.31/–1.68
(331)	+2.34/–1.17	+0.02/–0.01/+2.36/–1.18	+0.48/–0.44/+2.30/–1.25	+1.00/–0.48/+2.31/–1.68

Structure and Energetics

The energetic profile of H2 dissociation is summarized in Table , where all adsorption energies are referenced to the energy of H2(g) and CeO2 surface and were calculated according to

To the best of our knowledge, the activation energy of H2 dissociation was only reported for the most stable (111) termination. Our results for this surface, and in particular our activation energy (0.75 eV), are in perfect agreement with the ones reported by Fernández-Torre et al.[36] (0.76 eV), García-Melchor and López[34] (0.85 eV), Negreiros et al.[33] (0.72 eV), and Werner et al.[25] (0.78 eV). The activation energy increases in the series: (100)-O-(b) < (100)-Ce < (100)-O-(a) < (110) < (221) < (331) < (111) (Figure and Table ). These results show that the (111) termination exhibits the worst catalytic properties for H2 dissociation in the series and that it would be possible to reduce up to 5 times the energetic barrier by using (100)-O- or Ce-terminated surfaces instead of (111). For the (100) surface, it is important to note that the recent work of Capdevila-Cortada and López[65] has shown that the structural configuration of surface oxygen depends on the temperature. The (100)-O-(a) model is the most stable configuration and hence the most probable at low temperature (<100 K). However, for higher temperatures, the other configurations of surface oxygen can be accessible, and their probabilities increase with the temperature. Because these surfaces are less stable than the first one, their reactivity is expected to be higher, that is, with a smaller activation energy, as it was found in the (100)-O-(b) model.

Figure 3

Energetic profile (in eV) for the heterolytic dissociation of H2 molecule on CeO2 surfaces.

From an energetic point of view, the surfaces can be distinguished according to two groups: group I made of (100)-Ce, (100)-O-(a), and (100)-O-(b) terminations and group II regrouping (110) and (111) as well as step models (221) and (331). Group I and group II differ according to (i) the stabilization of the dissociation product (H+, H–) and (ii) the existence of an MS. Indeed, the surfaces of group I show a negative reaction energy (−0.20 to −0.64 eV), reflecting the fact that (H+, H–) is thermodynamically favored, whereas the surfaces of group II lead to an endothermic reaction (0.34–0.77 eV). The forward activation energy required to form the (H+, H–) intermediate lies between 0.16 and 0.35 eV for group I and between 0.53 and 0.75 eV for group II. Moreover, the backward activation energies for the surfaces of group II are quite small (between 0.01 and 0.11 eV), indicating the easiest reversibility of H2 dissociation on these surfaces as opposed to group I surfaces that have large backward energetic barriers because of a better stabilization of (H+, H–). Furthermore, only the surfaces of group II pass through an MS state during H2 dissociation. For the sake of completeness, the stability of the final OH/OH product is found to be −2.35 to −3.85 eV for group I and −2.18 to −3.19 eV for group II. This indicates a clear driving force for the formation of hydroxylated surfaces accompanied by a reduction of Ce4+ to Ce3+. The barrier for this process is reported for the (111) termination to be of 0.26 eV,[37] indicating a possible rapid, irreversible evolution of the hydride intermediate to the hydroxylated product, which is 2.84 eV more stable.

An analysis of the transition structure (TS) geometry shows similarities among the terminations (see Table and Figure that illustrate the case of the (110) termination). In TS, the H2 molecule is elongated, pointing to a Ce and an O in the surface. The distances H–H, Ce–H, and O–H are intermediate between the H2* and (H+, H–) intermediates: H–H distances span from 0.940 to 1.153 Å; Ce–H distances span from 2.434 to 2.647 Å, and O–H values span from 1.243 to 1.353 Å depending on their coordinations. Interestingly we found a relationship between the TS geometry and the related activation energy. More precisely, TS of low energy is associated to a geometry close to the H2* state (short H–H bonds and large Ce–H and O–H distances; Figure ) and correspond to group I, that is, the (100)-terminated models.

Table 3

Distance (in Å) of H–H, Ce–H, O–H, O–O, and Ce–O for TS Structure

	(100)-Ce	(100)-O-(a)	(100)-O-(b)	(110)	(111)	(221)	(331)
dH–H	1.072	0.960	0.940	1.103	1.153	1.106	1.126
dCe–H	2.510	2.498/2.434	2.647/2.580	2.335	2.242	2.310	2.306
dO–H	1.243	1.340	1.353	1.179	1.141	1.177	1.162
dO–O	2.687	3.842	2.969	2.678	4.154	3.006	2.963
dCe–O	2.218	2.384/2.248	2.228/2.301	2.688	2.742	2.646	2.604

Figure 4

Charge density difference of TS on the (110) surface. Blue and green isosurfaces show an electronic density gain and depletion, respectively.

Figure 5

(a) H–H, (b) O–H, and (c) Ce–H distances in TS structure as a function of H2 dissociation activation energy.

Electronic Structure

From an electronic point of view, the different steps involved in H2 dissociation were clearly identified in the density of states (DOS). All the structures show similar features, and we display (100)-O-(b) and (331) in Figure , with the projections on the Ce–H and O–H atoms for the sake of clarity. The projections on all the Ce, O, and H atoms are presented in the Supporting Information. First, H2 physisorption leads to the appearance of a narrow H2 molecular band in the valence region. In TS, the splitting of this H2 band into Hδ+ and Hδ− bands is observed. More precisely, the proton and hydride states are confined in two distinct regions, for all the surfaces studied. Indeed, although the position of the Hδ+ band is localized around −6 eV/EFermi, the Hδ− band is shifted to higher energies. More interestingly, only the surfaces of group I show an overlap between the TS and the slab bands. It has been shown that this overlap allows for a reduction in the transition structure energy in metals for the hydrogen evolution reaction.[66] Our results point in the same direction: the surfaces of group I show smaller activation energy and overlap between the levels of H2 in TS, whereas the surfaces of group II exhibit larger activation barriers, with no overlapping observed. For the product of dissociation (H+, H–), the hydride band becomes the highest occupied energy level. The main consequence is the increase of the Fermi level of the system. Also, the overlap between the oxygen (cerium) and H+ (H–) states indicates the presence of O–H+ (Ce–H–) bonds.

Figure 6

PDOS of the slab, H2*, TS, and hydride–proton pair (H+, H–) for the (100)-O-(b) (left) and (331) (right) surfaces. Only the cerium and oxygen bonds to hydrogen are projected. The PDOS of cerium, oxygen, hydrogen molecule, hydride (×3), and proton (×3) are depicted in blue, red, gray, fuchsia, and green, respectively. PDOS are given in arbitrary units, and the Fermi level is depicted with a dashed line.

Different behaviors between the group I and group II surfaces were observed. In particular, the surfaces of group I (with a small activation energy) show no or a slight increase in the Fermi energy during the dissociation (Figure ), contrary to the surfaces of group II (requiring a large energetic barrier) which involve a significant increase in the Fermi energy (Figure ). More details can be found in the Supporting Information, together with the projected density of states. Moreover, although the surface cerium atom shows very similar bands, oxygens can be differentiated in function of their structural configuration: twofold coordinated oxygens present a narrow peak around −5.5 eV/EFermi after H2 dissociation, whereas the peak is located near −6.5 eV/EFermi for threefold coordinated oxygen.

Bader charges analysis (Table ) indicates the electronic transfer during the dissociation process: first (H2* species), a preactivation of H2 molecule with a slight polarization of H–H bond is observed. In the TS species, the polarization of the H–H bond increases, leading to the formation of a tight ion pair, and finally, the system evolves to the formation of a hydride and proton, well-identified with their Bader charges: (H+, H–) species. Moreover, the oxygen involved in the formation of the hydroxyl group shows an important electron gain: approximately +0.5|e–| with respect to the same oxygen in the slab.

Effect of Temperature

The effect of temperature on the energetic profile was considered by calculating Gibbs free energies for T = 300 K (Table ). First, we note that temperature affects the stability of a species, and more precisely, the increase of the temperature goes along with the destabilization of all the adsorbed species. Indeed, at 300 K, a shift toward higher energy is observed: around +0.45 for TS; +0.50 for H2*, MS, and (H+, H–); and approximately +0.70 eV for (OH, OH). Therefore, the consideration of both enthalpy and entropy contributions shows that the heterolytic product (H+, H–) is only stabilized, that is, with a negative adsorption energy, on the (100)-Ce surface, whereas, without taking into account the temperature, this species was stable on the (100)-O-(a), (100)-O-(b), and (100)-Ce terminations. Although temperature induces very small changes (less than ±0.04 eV) in reaction energies, it reduces systematically forward and backward activation energies up to 0.09 eV. This decrease of the energetic barrier comes from the stronger destabilizing effect on H2* and (H+, H–) than on TS (Table ). Indeed, the imaginary frequency of the transition structure is not taken into account in the Gibbs free energy calculation, resulting in lower vibrational effects for the transition state. Moreover, at 300 K, for the (110), (111), (221), and (331) slabs characterized by a smooth potential energy surface between the TS and (H+, H–) states, the transition structure becomes more stable, that is, with a smaller adsorption energy, than the heterolytic product (H+, H–). This could indicate, as suggested by Negreiros et al.,[33] that H2 dissociation on ceria surfaces follows different pathways (from a reaction coordinate point of view) related to temperature.

IR Spectra

The vibrational frequencies and IR spectra of the H2 dissociation product (H+, H–) were computed for all the surfaces studied, except for (111) where the IR spectrum modeled corresponds to the MS because of the absence of the (H+, H–) state. The vibrational modes are distributed in three distinct regions (Figure ): O–H and Ce–H bending are characterized by low frequencies (between 450 and 900 cm–1), Ce–H stretching modes are comprised between 950 and 1250 cm–1, and O–H stretching shows higher frequencies with the wavenumbers in the range of 3100–3900 cm–1. The wavenumbers and their relative intensities are given in Table .

Figure 7

Computed IR spectra of the (H+, H–) species for (100)-Ce, (100)-O-(a), (100)-O-(b), (110), (221), and (331). For (111), the spectrum corresponds to MS because no (H+, H–) state was found for this surface. Intensities are given in arbitrary units.

Table 5

Wavenumbers (in cm–1) and Their Relative Intensities (in Brackets) of Ce–H– and O–H+ Stretching Modes of (H+, H–) Species for (100)-Ce, (100)-O-(a), (100)-O-(b), (110), (221), and (331) and MS Species for (111)

	(100)-Ce	(100)-O-(a)	(100)-O-(b)	(110)	(111)	(221)	(331)
	1095.3 (1.00)	1220.3 (1.00)	991.4 (1.00)	1172.4 (1.00)	933.1 (0.11)	1184.0 (1.00)	1184.7 (1.00)
νO–Hstretch	3715.1 (0.69)	3817.5 (0.83)	3151.4 (0.77)	3494.6 (0.85)	1769.3 (0.38)	3740.3 (0.24)	3710.7 (0.41)

For all the surfaces studied, the ν(Ce–H) vibration shows systematically the highest intensity. To the best of our knowledge, there is no experimental observation reported for the Ce–H vibrational mode. Previous works of Wu et al.[28] on CeO2 and CeH3 as well as a recent review of Copéret et al.[67] on surface hydrides suggest that Ce–H vibrations are expected in the IR spectral region in the range of 1500–2000 cm–1. According to our results, this hydride–cerium stretch region is found between 900 and 1300 cm–1. However, from an experimental point of view, the observation of this hydride vibrational mode is not straightforward, first because of the poor stability of the (H+, H−) product, and second because this is also the vibrational region of other adsorbed molecules as methoxide species (ν(C−O) stretching mode).[68,69] As regards the ν(O–H) vibrations, the wavenumbers are in good agreement with previous experimental[28,30,31,69,70] and theoretical[25,35] works, showing a peak in the region between 2700 and 3800 cm–1 as reported by Lustemberg et al.[71] The only exception is the (111) MS with a lower wavenumber (ν(O–H) = 1769.3 cm–1). Note that the presence of hydrogen bonds on the (100)-Ce and (110) surfaces reduces significantly the wavenumber of O–H stretching, which is in line with similar observations on the (111) surfaces previously reported by Lustemberg et al.[71] and Fernández-Torre et al.,[72] and explains the wide region observed for O–H stretching.

Descriptors of Reactivity

Capdevila-Cortada et al.[73,74] have defined some reactivity descriptors for ceria. For H2 dissociation, these descriptors were based on the basicity of oxygens [with the position of the O (2p) band] and on the reduction energy. We have explored the correlation between the structural and electronic parameters to derive descriptors that account for reactivity. We found that the activation energy displays a linear correlation with (i) the bond length between the cerium and oxygen surface sites in the bare relaxed slab and (ii) the difference in the Fermi energy between the dissociation product (H+, H–) and the bare slab. Indeed, a short Ce–O distance is associated with a small energetic barrier (Figure ) and is compatible with the TS structure with a low activation energy, as discussed in Section . Besides, a low activation energy is observed when the Fermi energy of the product (H+, H–) is close to the one of the slab (Figure ), which is consistent with an overlapping of the TS levels of H– with the slab, as happens in metals with the H evolution reaction.[66]

Figure 8

Activation energies as a function of Ce–O bond lengths for the bare slab.

Figure 9

Activation energies as a function of difference in Fermi energy between the dissociation product (H+, H–) and the bare slab.

Conclusions

In summary, we have shown that the topology of the surface plays a crucial role in H2 dissociation on ceria. In particular, the (100)-terminated surfaces reduce up to 5 times the activation energy, from 0.75 eV (111) to 0.16 eV (100). Moreover, these surfaces lead to the stabilization of the heterolytic product, forming pairs of hydride/proton species that are found to be thermodynamically stable in the (100)-terminated slabs. Harmonic IR vibrational spectra of the product were also computed showing the expected vibrational region of hydride species on ceria surfaces in the region of 900–1300 cm–1. Finally, we have shown that both the Ce–O bond length and Fermi energy can be used as the descriptors of reactivity for prediction purposes.

Methods and Models

Methods

All our calculations were based on the DFT framework combined with the projector augmented wave method,[75,76] as implemented in the Vienna ab initio simulation package[77−80] (VASP)—version 5.4.4. The core electrons were kept frozen and replaced by pseudopotentials (standard pseudopotentials for cerium and hydrogen and soft one for oxygen), whereas valence electrons (Ce: 4f15s25p65d16s2; O: 2s22p4; and H: 1s1) were expanded in a set of plane-wave basis functions with a kinetic energy cutoff of 300 eV. The generalized gradient approximation[81] approach was used for the exchange and correlation potential with the Perdew–Burke–Erzenhof (PBE) functional.[82,83] Dispersion forces were not considered because a previous work of Fernández-Torre et al.[36] has shown that they can be neglected for hydrogen in interaction with ceria. On-site Coulomb interaction corrections (DFT + U) following the Dudarev et al. approach[84,85] were employed for the cerium 4f electrons to partially correct for the self-interaction error. Ueff was set to 5.0 eV for the cerium 4f orbitals, in line with the previous works using typical U values in the range of 3.0–5.5 eV.[32,60,86−93] The Brillouin zone was sampled with the Monkhorst–Pack[94]k-point mesh with the k-points spaced 0.05 Å–1. The energies were converged to 10–4 eV in the self-consistent field and a conjugate-gradient algorithm was used to relax the atomic position until the forces acting on each atom were less than 0.01 eV/Å. Spin-polarized (unrestricted formalism) calculations were performed for open-shell systems. Bader charge analysis[95−98] was performed to analyze the electron and spin densities.[61−63]

The transition-structure (TS) structures were located by employing climbing image nudged elastic band (CI-NEB) algorithm.[99,100] A minimum of four images along the reaction pathway with a spring constant between the images set to 5.0 eV/Å2 and a convergence threshold on the forces acting on each ion of 0.02 eV/Å were used for all the CI-NEB calculations. Then, in a second step, the TS structures were refined by using the improved dimer method.[101,102]

Frequency calculations were performed (i) on each optimized structure to check that the related structures match to a minimum, (ii) and on each TS structure where only one imaginary frequency was obtained, demonstrating that the related TS structure correspond to a first-order saddle point. They were carried out with the finite displacement method as implemented in VASP. The Hessian matrix is determined from the harmonic force constants calculated as the numerical derivatives, with the atomic displacements set to ±0.015 Å. Only the top two atomic layers and the adsorbate were allowed to be displaced during the vibrational frequency run, with a threshold of 10–6 eV in the ionic loop. Harmonic frequencies and their associated normal modes were obtained by the diagonalization of the Hessian matrix, obtaining an accurate reproduction of low-energy vibrational modes (in the 500–1000 cm–1 range) and O–H stretching.[103]

The IR vibrational intensities, I, are proportional to the square of the first derivative of the dipole moment along the z-direction, dμ/dQ. Thus, the IR intensities were evaluated for each vibrational normal mode according towhere i is the ith eigenvector, is the mass-weighted coordinate matrix, μ is the dipole moment along z, and r is the atomic Cartesian coordinate.

Zero-point energy (ZPE) was calculated as followswhere ν is the frequency associated to the vibrational normal mode N.

Moreover, the effect of temperature was considered by calculating the Gibbs free energy: G = H – TS. If we assume that the expansion term pV in H = U + pV is negligible, which is a reasonable approximation in solid state,[104,105] then the Gibbs free energy can be approximated as the Helmholtz free energy

Therefore, for an adsorption system, a common and reasonable approximation is to consider only the vibrational contributions in the partition function,[104−106] leading to the following expressions for internal energy and entropy: U(T) = EDFT + EZPE + Uvib(T) and S(T) = Svib(T), respectively. More details about the thermostatistic equations used and the Gibbs free energy calculation can be found in the Supporting Information.

Slab Models

Ceria is a semiconductor material crystallizing in a cubic fluorite structure (Fm3m). The optimized cell parameter a calculated with PBE + U (U = 5 eV) was slightly overestimated by about 1.5% (aPBE+ = 5.495 Å) compared to the experimental value[107,108] (aexp = 5.411 Å) but is in good agreement with previous similar theoretical studies[60,89,91,109,110] and was used to build the slab models. The three stoichiometric low Miller index surfaces as well as step models were investigated: type I (110); type II (111), (221), and (331); and type III (100) according to Tasker classification.[111] Because the type III surfaces have a nonzero dipole moment normal to the surface, half of the surface atoms were moved to the bottom part to quench the dipole moment. The (100) surface can be O- or Ce-terminated. If it is well-known that the O-terminated surface is more stable than the Ce-terminated one,[64] recently Capdevila-Cortada and López[65] have reported the complexity of the (100) surface. In this context, we have built three models for this termination: one model for Ce termination—(100)-Ce; and two models for O termination—(100)-O-(a) and (100)-O-(b) (Figure ). Supercells were chosen to have surface models with approximately the same surface area (see Table ). The slab thickness was chosen to converge the surface energy for each termination within 5 × 10–3 J m–2 (see the Supporting Information for further details). The vacuum layer thickness was set to 15 Å for all the slabs to avoid interactions between images. In the first step, the bare slab built by cutting the CeO2 bulk was fully relaxed, and in the second step the adsorbate and the top-half slab were relaxed, keeping the bottom-half slab frozen. The main characteristics (geometry, atom coordination, and surface energy) of the seven slabs are shown in Table , and the models used are given in the Supporting Information.