A Deeper Insight on the Stability of Water-Induced Reconstruction of Anatase (001) Surface

: TiO 2 anatase (001) surface has been indicated for many years as a potential system for water dissociation and hydrogen production. Surface reconstruction periodicity of TiO 2 anatase (001) in water is revised on the basis of the new water induced reconstruction model that accounts for dissociative water adsorption in the ﬁrst monolayer and self-assembling of surface hydroxyls. The study has been performed in the context of ﬁrst principles total energy calculations on the basis of state of the art Density Functional Theory. Different surface periodical structures have been studied and compared in terms of residual surface stress and surface reactivity. While a preference seems to emerge for the (2 × 3) surface reconstruction, there are indications that this conﬁguration might not occur spontaneously in bulk water.


Introduction
Bulk TiO 2 and TiO 2 surfaces of the various crystal phases (anatase, rutile, and brookite) are being considered since many years for a variety of applications with special mention to environment-and energy-related applications, such as the photo-decomposition of organic pollutants, solar cells, and solar-hydrogen production [1][2][3]. Titanium oxides are also important in surgical implants where the implant surface is oxidized and in contact with the bio-moieties in the human body [4]. All the mentioned systems and problems require the knowledge of the surface-water interaction processes because they imply the usage of either nano-sized TiO 2 or TiO 2 surfaces in water. For this reason, the atomistic scale characterization and modeling of the TiO 2 -water interface is being studied intensively [5][6][7][8]. Among the crystal phases, the metastable anatase form is known to be photo-catalytically more active than the ground state rutile phase, and becomes the most stable phase at the nano-scale. This is why the anatase is increasingly considered. While the (101) is the most stable anatase surfaces, the minoritary (001) has been claimed to be more reactive in water [7,[9][10][11][12][13][14] and this is the reason why many articles have focussed on the (001) anatase surface in water. Moreover, recent results have reported on the successful synthesis of TiO 2 particles with minority facets [15][16][17]. The reactivity of the (001) surface with respect to the water adsorption and dissociation has been debated for long time with no conclusive statements yet. While the bare (unreconstructed) TiO 2 (001) surface has been predicted as a dissociative reactive surface for water absorption and, therefore, a good candidate to obtain photo-reactive water splitting for energetic purposes, the experimental results are rather controversial and do not confirm such predictions [18][19][20]. The (001) unreconstructed (dry) surface is made of rows of under-coordinated atoms along the [100] direction formed by sequences of five-fold titanium atoms (Ti 5c ) and two-fold oxygen atoms (O 2c ) with alternating long and short Ti 5c -O 2c bonds and an inward relaxation of the topmost layer with respect to the unrelaxed slab ( Figure 1). There is a general consensus that such a surface should be highly reactive for water dissociation [7,[9][10][11][12][13][14] but this is in marked contrast with the experiments that show, quite surprisingly, that this is not the case [18][19][20][21]. Moreover, in-situ experimental findings have shown that stable (001) anatase surface shows a (1 × 4) reconstruction character in vacuum. Atomistic models, based on popular first principles techniques such as the Density Functional Theory (DFT), have been constructed starting from the well known ADM model [22] that is able to reproduce the (1 × 4) surface reconstruction observed in vacuum. However there have been increasing difficulties to explain the experimentally observed surface intertness for water dissociation using the ADM model [21,23]. The reactivity loss has been related to the decreasing of the surface stress with respect to the unreconstructed surface [11][12][13] but, while the ADM model predicts a reduced surface stress due to the addition of TiO 2 molecules on the surface, this appears to be not sufficient because some residual surface reactivity, that is not confirmed by experiments, has been found [23]. Moreover it has been shown also that point defects, such as vacancies, oxygen ad-atoms or interstitials, might even enhance the ADM reactivity [23]. For this reason the ADM model has been further modified proposing that the observed surface reactivity might be related to the oxidized or reduced state of the surface [21]. Regarding the behavior of TiO 2 (001) surface in water, some authors have suggested that water patterned TiO 2 (001) surface might compete with the ADM models from the stability point of view [20,[24][25][26][27]. The present authors have recently revised the properties of water patterned (001) anatase surface in a large supercell [27] and two main new findings emerged: (1) Dissociative water adsorption forms a ridge-terrace geometry similar the one of the ADM model. This surface modification is accompanied by the formation of a stabilizing network of hydrogen bonds so that the water induced reconstruction (WIR) geometry shows an apparent (2 × 4) periodicity (see Figure 2); (2) such a water induced surface reconstruction reduces the surface stress at a larger amount than the ADM model and the surface reactivity disappears. However, the above mentioned WIR model was not studied and discussed in terms of the stability hierarchy of the possible various surface periodicities, even though previous studies have indicated the (1 × 3) reconstructed surface as the most stable one [24]. Then, in the present article, we try to go deeper into this aspect concerning the newly found stable ridge geometry of the WIR model.

Methods
All the calculations have been performed in the context of state of the art DFT with plane waves basis for both the electron density and wave-functions and periodic boundary conditions (PBCs). We have employed ultrasoft pseudopotentials [28] and the generalized gradient approximation (GGA) as implemented by the Perdew-Burke-Ernzerhof (PBE) scheme [29]. This computational approach has been already successfully employed in the previous literature revealing that PBE is the best GGA exchange-correlation scheme for the TiO 2 -water system [30]. Some technical aspects are worth to be mentioned such as the energy cutoff values employed for the wave functions and the electron density (60 Ry and 600 Ry respectively) and the threshold values employed for convergence of the atomistic optimization procedure, namely 10 −4 Ry/Bohr for residual forces and 10 −8 Ry for the total energy. It is important to evidence the necessity to add a dispersion correction term to the total energy in order to represent correctly the water-surface interaction. Generally speaking, this term is always expected to improve the energetics even though its inclusion causes a slight shortening of bond lengths and lattice constants [31,32] that, however, remain within the DFT-PBE error bar. In the particular case of the TiO 2 -water system, we have recently shown that the inclusion of the dispersion term is mandatory to get a reliable description of the adsorption and the surface energy values for the reconstructed surface. For all these reasons, we have included in the Hamiltonian a dispersion energy term by employing the well known Grimme-D2 correction [33]. Special mention deserves to be devoted to the well known self-interaction problem affecting standard DFT resulting in a wrong description of the electron localization of 3d orbitals. This inaccurate description affects the electronic properties of the system and should be taken into account especially for transition metals such as Ti. A common method to correct the self-interaction error is the so called DFT+U scheme where the orbital occupancy is forced to integer values by the inclusion of the Hubbard term in the Hamiltonian [20]. However we have shown that, concerning the present system, the inclusion of the Hubbard term does not affect the energetic properties we are interested in because these properties, depending on the difference between the energy values obtained from different configurations (see below), are rather insensitive to the inclusion of the Hubbard term [27] . For this reason the Hubbard term has not been included in our model.
Further details of the adopted computational scheme, together with the comparison of the bulk anatase properties obtained with the present scheme and other ones from the previous literature, can be found in a recent article [27].
The anatase (001) surface has been modeled through a slab geometry with the surface orthogonal to the z axis. The anatase bulk crystal was cleaved to expose the (001) surface and then relaxed. The slab contains four layers that provide an adequate geometry to describe in a correct way the relaxation pattern of the topmost layers in a real surface [23,25,32,34]. Moreover the supercell size along the z-axis has been chosen to have a ∼40 Å-thick vacuum region in order to avoid artifacts due to the PBCs.
An important choice is to keep fixed just the atoms belonging to bottom-most layer in order to reproduce the bulk atomic positions. Indeed, the less rigorous choice, adopted by many authors, of keeping fixed two bottom layers does not guarantee the correct relaxation patten in the slab [27].
Because the aim of the present article is to give a deeper insight to the water induced surface reconstruction periodicity and to the relevant water adsorption energy values, we have adopted different supercell sizes depending on the periodical ridge-terrace geometry investigated, namely (2 × 2), (2 × 3), (2 × 4), and (2 × 6) where the indexes refer respectively to the ridge direction [010] and to the orthogonal one [100]. For all the supercell sizes and surface reconstructions we have adopted a rigorous Monkhorst-Pack k-point sampling scheme (MP) [35] with denser and more demanding grids than the ones typically employed in the recent literature: our reference was to keep the grid at least as dense as the 2 × 2 × 1 MP grid for a (4 × 4) surface geometry. In other words we have employed a 4 × 4 × 1 MP grid for the (2 × 2) slab geometry, a 4 × 3 × 1 MP grid for a (2 × 3) surface, a 4 × 2 × 1 for (2 × 4) surface, and lastly 4 × 2 × 1 for (2 × 6). Lastly, the artificial electric field across the slab generated by periodic boundary conditions (PBC) along the z-axis is handled by using a self-consistent dipole correction scheme [36]. All the calculations have been performed with the Quantum ESPRESSO (QE) package [37].

Results
The formation energy of the unreconstructed surface γ is calculated according to the following equation: where γ f , E slab , E f slab are the surface formation energy of the unrelaxed slab and the total energy values respectively of the relaxed and the unrelaxed slab, n TiO 2 is the number of TiO 2 units in the slab, E TiO 2 is the energy of the TiO 2 unit calculated in the anatase bulk structure, and, finally, A is the surface area [27]. The unreconstructed surface energy is, of course, independent on the surface periodicity employed, i.e., on the n value (with n = 2, 3, 4, 6) of the (2xn) supercell. At the present level of the theory and using the above specified MP sampling scheme in the reciprocal space, we get γ = 1.105 J/m 2 , which is quite a large value reflecting the reactivity of the ideal unreconstructed surface.
Concerning the (2xn) water reconstructed surface, the surface energy can be calculated provided the average adsorption energy of dissociated waters is known: Indeed, if this is calculated as where E wet (2xn) and E dry (2xn) are the total energy values of respectively the water reconstructed surface and the unreconstructed one, n H 2 O is the number of water molecules dissociatively adsorbed and E H 2 O is the energy of an isolated water molecule, the surface energy can be obtained through: The resulting values of the average adsorption surface energy are reported in Table 1. The above values, which have been obtained in vacuum, can be further corrected by taking into account the dispersion energy contribution that should be considered in water. In bulk water, indeed, the adsorption process and energy, and the corresponding surface energy, must contain also the contribution from the dispersion energy change due to the water molecules that move from the bulk water to the surface where they are dissociatively adsorbed. This aspect has been detailed in a recent article [27]: According to this model, an estimate of such correction can be done with the following assumptions: • Provided that the total energy of the real system is obtained by summing the DFT and the dispersion energy terms, E = E DFT + E d , the interaction between the water molecules in bulk water is accounted for by the dispersion energy only (E DFT the dispersion energy change occurring in a real system when one water molecule is dissociatively adsorbed from the bulk water to the surface is obtained in a Born Oppenheimer Molecular Dynamics (BOMD) run as detailed recently [27]; this contribution is supposed to be due only to the change of the dispersion energy interaction between water molecules and between the TiO 2 slab and water molecules (the dispersion energy change due to the TiO 2 re-arrangement is neglected); BOMD simulations were performed on on the anatase unreconstructed and fully hydrated (001) surface with the mixed Gaussian and plane wave (GPW) method as implemented in the CP2K package [38]. The details and technicalities of the simulations can be found in a recent article that was mostly focussed on the study of the water structure [39]; • the contribution to the dispersion energy change due to the interaction between the adsorbed (and dissociated) water and the TiO 2 slab is measured from the ground states of the water reconstructed and the unreconstructed surfaces; • the dispersion energy of one water molecule in bulk water is measured and averaged from BOMD of bulk water at room temperature; for this term we have obtained E d H 2 O ≈ −0.0245 ± 2 * 10 −3 Ry; • the contribution to the adsorption energy coming from the interaction between the slab and the bulk water above it is dominated by the medium-long range interactions and thus by the dispersion part of the energy.
Using the above assumptions it can be shown [27] that the adsorption energy, appropriately corrected with the dispersion terms that should be present in bulk water, casts into: where (H 2 O) m(d) stands for molecular (dissociated) water molecules and E d are respectively the dispersion energy of the TiO 2 slab interacting with q water molecules in bulk water and the dispersion energy term of the two adsorbed and dissociated water molecules interacting with p water molecules in bulk water. On the basis of this simple model, we have calculated the dispersion corrected values of the average adsorption and surface energy for (2xn) reconstructed surfaces with n = 2-4 as listed in Table 2: The surface energy values corrected with the dispersion energy contribution from the bulk water show a new stability hierarchy where the second most stable surface reconstruction is the (2 × 4) one instead of the (2 × 2) one. Moreover, while the surface energy values increase with respect to the one in vacuum, the difference between (2 × 3) and (2 × 4) is reduced. At this point, we have employed the (2 × 6) surface reconstruction in order to reproduce in a doubled supercell the (2 × 3) surface reconstruction. We have attempted to obtain a spontaneous dissociative adsorption in the middle of the terrace so as to reproduce a (2 × 3) periodicity and we have checked several different initial configurations. However we have only obtained molecular water adsorption on the terrace, as reported in Figure 3.

Discussion
Two main results of the previous section are worth to be analyzed and discussed: first the stability hierarchy between the three WIR surface geometries that have been considered, namely (2 × 2), (2 × 3), and (2 × 4), is changed when the missing part of the dispersion energy, the one that is typically due to the bulk water above the surface, is included. Hence, the results obtained in vacuum might change in bulk water as evidenced by the change of the positions in the stability hierarchy between (2 × 2) and (2 × 4): In bulk water, indeed, the (2 × 4) surface is more stable than the (2 × 2) that, on the contrary, appears to be more stable in vacuum. This behavior is related to the different values measured for the ∆E d slab−bulk water energy data shown in Table 2: the dispersion energy change, due to the interaction with the bulk water, between the reconstructed and the unreconstructed surfaces increases with the surface area (also this term per unit area increases with the area as well). Then, according to the Equation (4), the dispersion energy of bulk water E d H 2 O contributes to lower the adsorption energy. The surface energy is obtained from the corrected values of the average adsorption energy by dividing for the surface area and this reduces the positive correction term (that lowers the absolute value of the adsorption energy) at larger extent in larger areas. In other words, the correction per unit area is greater for (2 × 2) than for (2 × 4) thus resulting in the inversion of stability. It is interesting to note that the correction due to the dispersion energy results in closer surface energy values between (2 × 3) and (2 × 4). Moreover, the evidence that the adsorption in the (2 × 6) reconstructed surface is molecular also in the middle of the terrace, evidence that we have got after various attempts to obtain a dissociative adsorption, indicates that the (2 × 3) reconstruction should be, most probably, thermally activated because we expects the existence of an energy barrier for a further dissociative adsorption along [100] as soon as the ridge is already formed.
Lastly, because the adsorption on the (2 × 3) and (2 × 4) water reconstructed surfaces is molecular [27], it is plausible that in bulk water the hierarchy could be further changed if one takes into account the molecular adsorption energy occurring on the available sites. Indeed also molecular adsorption events might release some energy per unit area due to the small, but not negligible, surface relaxation. For this reason we have sampled all the possible adsorption sites on the (2 × 3) and (2 × 4) surface to check whether such inversion could occur in bulk water.
The adsorption sites sampled are reported in Figure 4 and the molecular adsorption values obtained in vacuum are reported in Table 3.   Table 3 are larger for (2 × 4) than for (2 × 3) and this is consistent with the lower residual stress found in (2 × 3). This circumstance, anyway, tells us that the surface relaxation induced by molecular adsorption onto (2 × 4) should be larger than the one occurring on (2 × 3). In addition, considering the (2 × 2) case, we see that although this surface evidences the lowest residual tensile stress, the surface energy is the largest one occurring in bulk water teaching us that the further inclusion of the dispersion energy correction makes the relationship between the residual surface stress and the surface energy in bulk water not that straightforward as in vacuum.
In addition, the inclusion of the dispersion energy correction would presumably play the same role also for molecular adsorption as the one found for dissociative adsorption. Considering the close values of the surface energy in bulk water for (2 × 3) and (2 × 4) water induced reconstructed surfaces, this circumstances, together with the presumable existence of an activation energy for dissociative water on a (2 × 6) surface, makes the question of finding the most stable WIR configuration in bulk water still open.

Conclusions
In summary we have revised the stability hierarchy of water induced reconstruction periodicities of anatase (001) surfaces on the basis of the newly found (2xn) geometry. Using the dispersion energy correction, we have shown that in bulk water the stability hierarchy found in vacuum is changed with the (2 × 4) surface reconstruction being more stable than the (2 × 2) one. The (2 × 3) surface appears to be the most stable one also in bulk water but the surface energy difference between (2 × 3) and (2 × 4) decreases with respect to the vacuum case. The molecular adsorption energy values on both the terrace and the ridge sites are larger (in absolute values) on the (2 × 4) than on the (2 × 3) in accordance with the surface energy values and the residual surface stress. This means that the molecular adsorption induces a larger relaxation pattern onto the (2 × 4) than onto the (2 × 3) surfaces. Thus we can argue that, if one considers the dispersion energy involved in the molecular adsorption, the same surface energy correction behavior could occur for molecular adsorption (to a lesser extent perhaps) as the one found for dissociative adsorption. Another point is that the water adsorption onto the (2 × 6) surface is not spontaneously dissociative, meaning that some activation energy must exist. For all these reasons, we can state that the real stability hierarchy of WIR anatase (001) surface in bulk water should not be considered as a closed question yet. Free energy measurement in the context of ab-initio molecular dynamics simulations could help to clarify this important issue.