Investigating Surface Properties and Lithium Diffusion in Brookite-TiO2

Surfaces properties of TiO2 in the brookite phase and the lithium diffusion are studied using density functional theory (DFT) and interatomic potential simulations. Simulations predict that the brookite surfaces have a higher intrinsic Lewis acidity compared to the other polymorphs due to the presence of four coordinated Ti atoms on the surface in contrast to the most stable surfaces of anatase and rutile which have five coordinated Ti surface atoms. The surface reactivity of the brookite is then expected to be higher for catalysis. The (210) surface is the most stable calculated surface and some experimental morphologies have been revisited. Regarding Li intercalation, the migration occurs along the c-direction and open channel surfaces are desired, therefore particles exposing the (001) surface. Those differences on chemical/physical properties highlight the importance to control the nanoparticle shapes for a given application.


Introduction
The anatase and rutile phases have been widely studied 1 because they are often used in photocatalysis, 2 as a catalyst support, 3 in dye-sensitized solar cell 4 and for medical application (the photocatalyst can be seen as the bactericidal, 5 the oxide surface of Ti is used for biocompatibility). 6 There are fewer experimental studies on brookite despite its potential for use in photocatalysis, [7][8][9] catalysis, 10,11 for photovoltaic applications [12][13][14] and lithium rechargeable batteries. 15,16 Brookite is a metastable phase of stoichiometric TiO 2 and is less common than the other phases, anatase and rutile. Experimental studies have derived the thermodynamic properties from 298 K 17 to high temperatures. 18 They confirm that rutile is the most stable phase followed by brookite and then by anatase.
Although studies have succeeded in synthesizing brookite, most of them also obtain the other phases, anatase and rutile. Furthermore, as the occurrence of the exposed surfaces is highly dependent on the synthesis conditions under which brookite is prepared (e.g., hydrothermal or hydrolysis), it is difficult to assess the relative stabilities of the different surfaces from experiment. This is further exemplified by very different morphologies that have been observed such as tablet/sheet, [19][20][21] stick, 9,22 pseudo-cubic 23 and spherical. 8 Some authors 19,20,24 have been able to identify the exposed surfaces using TEM (transmission electron microscopy) and XRD (X-ray diffraction). Moret et al. 25 have used X-ray methods and found that the dominant surfaces are (210), (111), (211) and (102). In contrast, Kominami et al. 24 have identified the (121) surface along with the (120) or (111) surface. Shibata et al. 20 have analysed brookite with a plate-like habit in which the major surface was (010). Finally, Pottier et al. 19 have obtained nanocrystals of brookite dominated by a basal face, supposed to be an unusual (301) facets with only two exposed surfaces surrounded by four (111) surfaces.
To better link the structures of brookite surfaces with their stabilities and reactivities, we aim to generate reliable models to express the shape of the particles regarding the experimentally observed surfaces. We combine DFT (density functional theory) and force field methods. The diffusion path of lithium in brookite using DFT methods is compared to the results by Arrouvel et al. 15 using force field methods. We focus with DFT on the (210), (100) and (010) surfaces but also on the (301) and (111) surfaces that have been identified by Pottier et al. 19 A comparison of brookite with previous simulations on TiO 2 polymorphs [26][27][28][29][30][31][32][33][34][35][36][37][38][39][40] is also addressed to enable us to infer the potential of brookite for catalysis and energy materials.

DFT calculations
Total energy calculations are performed within the density functional theory (DFT) implemented in the Vienna ab initio simulation package (VASP). [41][42][43] Different functionals have been tested: the generalized gradient approximation (GGA) of Perdew and Wang (PW91), 44 the revised Perdew-Burke-Ernzerhof functional (RPBE) 45 and local density approximation (LDA). 46 PW91 is the most accurate on the bulk study (see Results section) and is mainly used for this study. Another approach is to use DFT+U 47 which combines the GGA functional with U, the site Coulomb parameter, for the electronic correlation for d orbitals. U is a value adjusted by hand to reproduce some experimental properties. The influence of the U value on some properties is tested. To solve the Kohn-Sham equations, VASP performs an iterative diagonalization of the Kohn-Sham Hamiltonian via unconstrained band-by-band minimization of the norm of the residual vector to each eigenstate and via optimized charge density mixing routines. The convergence criterion for the electronic self-consistent cycle is fixed at 0.1 meV per cell. The eigenstates of the electron wave functions are expanded on a plane-waves basis set using pseudopotentials to describe the electron-ion interactions within the projector augmented waves (PAW) approach. 48 For total energy calculations with PW91, we use a cutoff energy of 312.5 eV. The optimization of the atomic geometry at 0 K is performed by determining the exact Hellman-Feynman forces acting on the ions for each optimization step and by using a conjugate gradient algorithm. A full relaxation of all atomic positions in the cell is performed until the geometric convergence criterion on the energy (1 meV per cell) is reached. A convergence on the k-point meshes using Monkhorst-Pack scheme 49 is also ensured on bulk structural optimization and it is even overvalued for the band gap calculations (see the values in the Results section).

ZPE contributions
The zero point energy (ZPE) contribution is calculated using the well-known lattice dynamics approach of Born and Huang, 50 where the vibrational frequencies are calculated by diagonalizing the dynamic matrix at different points in k-space and then integrated over k-space to yield the energy contribution. The problem is that the direct calculation of the vibrational frequencies as a function of k-space is CPU (central processing unit) time consuming when using ab initio approaches. However, model potentials 51 can achieve virtually the same level reliability, certainly when integrated over k-space, and have the advantage of being very quick and efficient. We used PARAPOCS program, 52,53 which was developed to perform free energy minimization of solids using numerical derivatives, as a consequence gives the vibrational frequencies and the different components of the free energy suitably evaluated over the whole Brillouin zone. This code has also been previously used in conjunction with DFT calculations modelling two low index surfaces of ZnO. 54 Interatomic potential methods The simulations are upon a Born model of solid where ions interact via long-range electrostatic forces and short range-forces including both the repulsions and the van der Waals attractions between neighboring electron charge clouds. 50 We have used interatomic potentials derived by Matsui and Akoagi 55 which were successfully used to model the surfaces of rutile 51 and adapted to model the rutilewater interface. 56 Those potentials have been described by Oliver et al. 51 and our atomistic simulations have been undertaken with the same code, METADISE (minimum energy techniques applied to dislocation, interface, and surface energies). 57 The potentials follow the relation: with q i , q j the charges of each ion, r ij the separation of the ion centers, and A, ρ, C the ion-ion parameters in the Buckingham relation. And they are given in Table 1. Those potentials do not have a shell model and their validity for brookite-TiO 2 has been compared to DFT on surface energy calculations.

Surface energy calculations
The surface energy in vacuum, , for a given crystal orientation (hkl), is given as follows: (2) G hkl is the Gibbs free energy of the (hkl) surface and the G bulk is the Gibbs free energy of the bulk normalized to the number N of atoms used in the supercell. A hkl is the surface area. We note that a slab of the supercell contains two (hkl) surfaces.
The calculations are performed at 0 K and neglect the zero point energy which is valuable considering that the ZPE difference between the surface and the bulk is negligible. Finally, the slab energies assume constant surface area. Therefore, we can write: G = E, with E the energy of the bulk or the slab. The equation 2 becomes: From the k-point grid optimization of the bulk energies (in Table 2) we can evaluate the k-point grid necessary for the calculations of each surface. To avoid the interaction between slabs along the z axis, we make a vacuum gap higher than 26 Å for each surface and ensure that the slabs have mirror symmetry (vacuum and slab sizes are summarized in Table 3). In each simulation, all of the atoms are allowed to relax. Another point to take into account is the eventual oscillation of the surface energy with the parity of the number of layers. When a surface has an oscillating surface energy, the minimum is observed for even numbers of layers and a maximum for odd numbers of layers, whatever functional or pseudopotential is chosen (for example on the (110) surface of rutile, 31,34,37,39 (110) and (100) surfaces of anatase). 58 In order to ensure that we had a reliable value of surface energy, we chose a relatively high and even number of layers in the case of the brookite and rutile phases, as the case of the anatase phase having been carefully detailed in a previous study. 26 The (101) and (001) surfaces of anatase do not oscillate and from a previous study 26 the required slab depth is known to have a reliable value for the surface energy. The details for each calculated surface are given in the following Results section. The value of the surface energy may vary with both slab thickness and with different functionals. However, the variation of energy is the same so it is important to keep the methodology consistent when comparing different surfaces. For that consistency, we choose the method from Arrouvel et al. 26

Gibbs-Curie-Wulff construction
The Gibbs-Curie-Wulff law 59 provides an approach determining the equilibrium morphology of an ideal crystal according to the following relationship, (4) where Γ hkl is the surface energy for the (hkl) orientation, d hkl is for the normal distance from the surface to the center of mass of the solid, and α is a constant independent from h, k, l.
The METADISE 57 code allows us to express the equilibrium shape of the oxide using the surface energies from DFT and atomistic calculations and indeed to compare it with predictions using simplified pair potentials. Previous theoretical studies 58,60 have shown that the shape at the thermodynamic equilibrium depends strongly on the working conditions. However, the first step is to consider the material under vacuum.

TiO 2 bulk
The crystal structure of brookite is orthorhombic (a = 9.184 Å, b = 5.447 Å, c = 5.145 Å). 61 Before constructing the surfaces of brookite, it is necessary to optimize the bulk cell from X-ray data 62 using the simulation approach described above. A comparison with the bulk crystal structure of anatase 63,64 and rutile 65 obtained from X-ray data are summarized in Table 2. From the convergence criteria on bulk calculations, the convergence on ionic forces is ca. 1 meV Å -1 .
The rutile phase is the densest of the three phases with the highest symmetry. The bulk cell ( Figure 1 From the calculated total energy E T , the relative stability of each phase can be deduced and compared with the experimental enthalpy differences. Effectively, the calculated total energy corresponds to the internal energy and can be simply linked to the enthalpy by neglecting the effect of temperature and zero point energy (ZPE). The experiment for the enthalpy data is given at ambient temperature and the extrapolation at 0 K does not reverse the order of stability between anatase and rutile. The calorimetric work has also shown that the standard entropy of the different phases are similar and hence entropic differences will not affect the order of stability (values are available for rutile and anatase). 17,67 The calculations show that anatase is more stable than rutile by 8 kJ mol -1 (∆E T per TiO 2 unit) and less surprising is more stable than brookite by 3 kJ mol -1 . We found that the ZPE, calculated empirically with PARAPOCS, 52,53 does not affect the order of stability. We have also tested the RPBE (revised Perdew-Burke-Ernzerhof) functional with a PAW pseudopotential and a cutoff of 353.6 eV; and still found that the anatase phase was the most stable one at 0 K. In addition, the cell parameters are less well represented compared to PW91. Using DFT+U 47 we found, in agreement with Vu et al., 68 that the rutile phase becomes the most stable phase with the increase of the U value. However, the choice of the Coulomb site U and the exchange site J on Ti is somewhat empirical. We chose to fit to the band gap of rutile phase at 3.0 eV using the values of U = 9 eV and J = 0.3 eV, but the fitted values give less reliable band gaps for the other phases, the anatase at 3.4 eV and the brookite 3.3 eV, i.e., the band gap of anatase is overestimated. A high value of U (U = 10 eV) has been reported with LDA by Person and Ferreira da Silva 69 in order to reproduce the band gap of rutile. From our results, U is not reliable enough to compare different polymorphs. Tables S2 and S3 in Supplementary Information (SI) section resume the calculated values of band gaps and cell volumes regarding different U values.
There are some theoretical studies on the bulk electronic properties but few reproduce the experimental band gap due to the choice of the functional. The LDA functional   underestimates the band gap (anatase: 2.0 eV; 70,71 rutile 1.5-2.0 eV), 71,72 and Hartree-Fock calculations overestimate it (anatase: 10 eV). 73 The result with a hybrid method B3LYP (Becke 3 Lee-Yang-Parr, functional combining three thermodynamic experimental parameters) 74 shows a better agreement (anatase: 3.7 eV). 27 The calculation close to experimental value is the one using LMTO (linear muffin-tin orbital) method on rutile (3.06 eV) 75 but complementary results are not available for anatase and brookite.
There are many experimental references in the literature to the band gap of bulk rutile, 3.0 eV 2 and anatase, 3.2 eV. 76 From our DFT calculations with GGA/PW91 functional, the band gap is also underestimated. The relative order is however consistent with experimental data. The band gap of anatase is 0.3 eV higher than the rutile. There is no experimental data on large crystals of brookite due to the lack of experimental pure specimens and that most experimental syntheses generate nanoparticles. Indeed, some studies of the band gap on nanoparticles of brookite by Koelsch et al. 13 give 3.4 or 3.26 eV, 20 which is still in agreement with the calculated relative order. Koelsch et al. 77 showed also that the band gap varies with the nanoparticle size. The authors do not know if the band gap variation is due to the direct or indirect transition. The direct transition increases when the size decreases. For the indirect transition, the variation is less pronounced but the tendency is similar on anatase and brookite, and more variable on rutile. Other authors report for nanoparticles of brookite band gap at 3.11 eV, 78 3.13 eV. 22 The study from Koelsch et al. 77 has suggested that the morphology plays a key role on the electronic properties and then, it becomes important to compare the different surfaces of each phase. Furthermore, as we want to explore whether brookite is more reactive and hence potentially a technologically more important material, we needed the simulations to be consistent with previous studies. 26 As the GGA/PW91 functional has been shown to be able to describe the surface structures and their reactivity we report the surface energies with this functional and PAW pseudopotentials. Some tests have been undertaken with DFT+U and the surface energy varies with the U value but the structure of the surface remains the same (not shown in the present paper).

Li diffusion in brookite-TiO 2
From a previous study, 15 we have identified the lithium pathway using force field. The lithium is diffusing along the c direction and two channels are interconnected. The optimized path is reproduced with 15 positions following the c direction in one of the channels and the connection between two channels, see Figure 2. For each calculation, the spin-polarization is considered, all ionic positions are relaxed and only the z coordinated is keeping fixed for Li to be able to evaluate the path and energies along the c direction. Some key Li positions have been fully relaxed in order to reach the fully relaxed stable intercalated lithium. We report an activation energy of 0.73 eV, slightly lower than the activation energy obtained with force field (E ac = 1.07 eV). 15 This difference is consistent with calculated activation energy barriers on TiO 2 −B, comparing the same methods (DFT and force field). 15,79 The most stable site at the center the unit cell relaxed at fractional coordinates 0.59 0.62 0.55 (see Figure 2b). This result agrees with the previous study from Arrouvel et al., 15 this site having been identified as the C site, with the Wyckoff notation 8c. The center of the unit cell is in fact a barrier between two c-channels, agreeing with force field methods, site M, Wyckoff notation 4b. The 6 oxygen atoms in the first coordination sphere are at distances 1.83, 1.86, 2.01, 2.03, 2.16 and 2.66 Å. We note that there are four equivalent positions. Another stable position is found at The DFT results are in excellent agreement with force field methods for geometrical properties. No magnetic moment is observed in the stable sites. During the diffusion, the total magnetic moment of the unit cell can reach 0.26 µ B in some positions close to the stable sites.

Brookite-TiO 2 surfaces
We will characterize the main surfaces which have been observed experimentally on particles from different methods. There is special interest in (301) and (111) surfaces because the synthesis of the brookite at low pH (pH = 1) from TiCl 4 induces the formation of well-defined nanoparticles exhibiting mainly those surfaces with the (301) identified as the basal face 19,77 by TEM. However, other studies from Ohno et al. 23 and Lin et al. 21 also demonstrate morphologies exposing the (210) surfaces. Resultantly, we describe other important surfaces, such as (210), (100) and (010) surfaces of brookite (Table 3). Brookite has the least symmetrical crystal structure than anatase and rutile, therefore more atoms are required in the simulation cells (higher number of layers and higher surface area) in order to have a good description of the surfaces. All the surfaces relax considerably, which means that a larger number of layers is required. Only the most stable terminations are presented. In those slabs, we can see that the surface energies before the ionic relaxation are twice the surface energies after relaxation. The convergence on ionic forces is ca. 10 meV on each atom for the surfaces.
We denote the acid Lewis sites by Ti V the five coordinated Ti atoms and by Ti IV the four coordinated Ti atoms. The basic sites generated are bridging oxygen µ 2 O.
(301) surface Modelling the (301) surface has been our first focus because it is the dominant surface of nanoparticles synthesized using TiCl 4 as discussed above. This surface has a large energy (after relaxation: 1.03 J m -2 ) in comparison to surfaces with lower Miller indices ( Figure 3 and Table 3). Before relaxation, the surface energy is  Ti IV and Ti V sites. In our surface unit cell (surface area of 99.26 Å 2 ), the four Ti IV atoms and the two Ti V atoms are non-equivalent (Figure 3) whereas on the (110) surface of anatase, all the Ti IV are equivalent.
We denote the number of broken Ti−O bonds which corresponds in fact to double the number of µ 2 O, being at 26.8 µmol m -2 . Effectively, Hu et al. 80 have shown on the various clays that the point zero of charge (PZC, also assimilated in that case to the isoelectric point) is linked to the number of broken bonds. Considering this hypothesis, we should be able to predict a relative order of the acidity of each surface.

(111) surface
The (111) surface of brookite has been characterized by different authors. 19,24,25 The slab unit cell is monoclinic (α = β = 90°, γ = 70.38°) and is more stable than the (301) surface (0.78 J m -2 ). The number of broken bonds is similar to the (301) surface, at 22.0 µmol m -2 . The surface has also different sites (Figure 4). In each surface unit cell, there are four surface Ti sites, one Ti IV and three different Ti V .
The presence of Ti IV indicates that this surface should also be reactive to small molecules. Effectively, a Ti IV should be able to adsorb two basic groups (such as OH or SH groups) to complete its bulk coordination of 6. Those sites are clearly more acidic than Ti V sites and as with the (110) surface of anatase, H 2 O or H 2 S molecules will dissociate on Ti IV sites. 60 Before relaxation, the surface energy is 1.74 J m -2 . This surface shows the six different Ti−O distances found in the bulk. The relaxation results in most of the distances being shortened by 0.1 Å, except for the two longest distances on the Ti V, which become even longer. The distances corresponding to the 2.081 Å bond lengths stretch to 2.211 Å. Thus the Ti V atoms relax so that they are more similar to Ti IV which will increase the Lewis site quality.

(210) surface
The (210) surface has also been characterized using electron diffraction. 20 The surface is the most stable surface ( = 0.71 J m -2 , in agreement with other theoretical work), 40,81,82 and has four Ti V sites in one unit   The (100) surface is terminated by Ti V and µ 2 O at a high density (11.68 µmol m -2 , Table 3) and is scarcely observed. 83 The number of broken bonds is then at 23.4 µmol m -2 .
The surface is more stable (at 0.95 J m -2 ) than the (310) and (111) surfaces and it is isostructural with the (112) surface of anatase. Figure 6 shows that each octahedral unit in anatase is aligned in the same direction whereas in brookite, the orientation of the octahedra changes every two layers. The Ti−Ti and Ti−O distances are regularly spaced. The Ti−Ti distance is 3.095 Å before relaxation and shortens to 2.800 Å after relaxation. All the Ti V atoms are identical. The µ 2 O are also similar with initial Ti−O distances of 1.927 and 2.081 Å which become, respectively, 1.799 and 1.884 Å after relaxation.
The (112) surface of anatase is also rarely seen. 84,85 Gao and Elder 85 observed the {112} face on anatase (and {103} face) synthesized at high pH. Furthermore, because the structure of the (100) surface of brookite is epitaxially matched to the (112) surface of anatase, it is not surprising that brookite has been observed at the surface of anatase. 83,86 If those surfaces are exposed, they should exhibit the same properties, but, as yet, there is no experimental data to confirm this fact.

(010) surface
The (010) surface (see Figure 7) has been detected by Shibata et al. 20 as the dominant surface in their sample. This surface is energetically the most stable (0.83 J m -2 ). This surface gets the highest number of broken bonds (27.7 µmol m -2 ). The acid sites are the Ti IV atoms with a high density (6.94 µmol m -2 ). That explains why the brookite from their study 20 is predominant and very hydrophilic especially when it has been stored under high vacuum or been exposed to UV light radiation. These conditions clean the surface of adsorbed molecules and the Lewis sites are free to accept water molecules. The Ti−O distances are 2.025 and 1.927 Å before relaxation and shorten at 1.829 and 1.858 Å after relaxation.
This surface has also been calculated by other authors. 66,87 They also found that the (010) surface has a lower surface energy than (100) surface.  (100) surface is not appearing, which is not surprising because it is not observed on the pure brookite sample. The main exposed surfaces are the (111) and the (210) (Figure 8). Those surfaces have been often observed. 19,20,24,25 The (010) is also present to a small extent, 5%. This surface has been strongly stabilized during the synthesis by Shibata et al. 20 as a flat particle exposing the (010) surface with the (210) surface at the side. The (301) surface appears with a small surface area (1%).

Gibbs-Curie-Wulff morphologies
Atomistic calculations using Matsui potentials give the almost same trend on the surface energy. We are able to simulate a higher number of surfaces (21 surfaces have been calculated, see Figure S1 in SI section). The surface    (101) and (212) surfaces, others can be even negligible such as the (100) and (011) surfaces being at less than 1%.
From Pottier et al. 19 study, the authors attribute the (301) surface as the basal from TEM analysis and the four sides surfaces belonging to (111) family. From our calculations, these surfaces are not the most stable ones. In addition, to obtain the tablet obtained experimentally, only two faces of (301) orientation should appear instead of four facets. To reproduce this morphology, the (301) surface needs to be strongly stabilized (i.e., the surface energy becoming the lowest). To obtain the same morphology as in the experiment (Figure 9), we have fixed a surface energy ratio (111)/(301) at 1.5. We note also that the edges are perpendicular to the basal surfaces, which correspond to the , and , planes perpendicular to the {301} and . However, another possibility more likely to explain the tablet shape is to consider that the basal plane is the (001) surface and the side surfaces are the (210) facets. A morphology with an angle of 80 degrees between the side facets and perpendicular to the basal plane has been observed by Ohno et al. 23 in some particles. This angle correspond the intersection of two (210) planes, measured at 80.6 degrees from our DFT calculations. The morphology obtained by Pottier et al. 19 is then revisited in Figure 9.

Surface electronic calculations of rutile, anatase and brookite phases
Electronic calculations require a high number of k-point meshes to converge the band gap to within 0.03 eV. Our study on the bulk of TiO 2 shows that it is important to keep the same methodology to compare the different phases and surfaces. We have reported in Table 3 the five surfaces of brookite. In addition, tests on the influence of the slab thickness on the band gap show that it modifies by 0.1 eV (the difference observed between the slab of (110) surface of rutile at 6 and 8 layers). It has been not possible to compare band gap of similar slab thickness due to steric effect and to reduce the (100) surface to 10 Å will not be appropriate (the slab is not thick enough and the surface energy is not converged).
The test on the rutile (110) agrees with previous experimental and theoretical studies 1,37,88 where the band gap of the (110) surface is similar to the bulk values. We found, for the slab with 6 layers of TiO 2 , a band gap at 2.06 eV (with k-point meshes at (15 × 10 × 1)) and at 1.97 eV for the 8 layers model which converges to the bulk value at 1.91 eV. Another study using hybrid DFT methods with a slab with 5 layers has reported that the gap is slightly smaller than the bulk by 0.2 eV. 89 The difference might come from the functional used or might be due to the odd number of layers. But by increasing the number of layers to 9, the band gap converges to the bulk value.
For the anatase phase, the band gaps of the (101) surface and the (001) surface are very different. The band gap, E g , on the (101) surface is higher than the bulk by 0.26 eV (E g = 2.48 eV) whereas the (001) surface is more conductive, being lower by 0.69 eV (E g = 1.53 eV).  The values for each surface of brookite phase are reported in Table 3. Only the (301) surface has a band gap slightly lower than the bulk value by 0.05 eV (E g = 2.27 eV). However, the (100) surface has a band gap close to the bulk, at 2.33 eV, which is similar to the relative order found by Beltrán et al. 66 The order for E g in Table 3 is given as follows: (301) < bulk ca. (100) < (210) < (111) < (010).

Discussion
Rutile is the densest phase of the three considered. ∆ f H° values 17 indicate that rutile is more stable at 298 K by 2 kJ mol -1 than brookite and by 5 kJ mol -1 than anatase. Brookite should therefore be more stable than anatase. ∆ f S° values are similar for rutile and anatase (ca. 50 J mol -1 K -1 ) 17,67 which will not influence the stability order (no value is available for brookite). In contrast with experimental measurements, we find from calculations and in agreement with various theoretical studies, 29,73,90 that the anatase phase is more stable than the rutile by 8 kJ mol -1 and the brookite phase more stable by 4 kJ mol -1 . The difference is small and the functionals (GGA PW91 and RPBE) may be at the origin from this discrepancy. The work by Muscat at al. 90 also shows that the choice of the basis set can influence the order, with a good basis set the anatase is still the most stable phase. Fahmi et al. 73 also found some differences when using HF, the anatase, being more stable, becomes less stable when the parameters are calculated with correlation density functional. The introduction of the ZPE correction does not change the stability order but emphasizes that we should be careful when comparing the energetics of different bulk polymorphs with DFT. However, we find that fitting the band gap on rutile at 3.0 eV with DFT+U, the rutile is the most stable phase but the corrections by U and J need to be the same to provide a meaningful comparison between the different phases. The band gap of anatase deviates from the experimental value when U and J are fitted to the rutile phase band gap. This is not surprising because the Ti environment is different for each phase. Titanium atoms are in an octahedral environment with a D 4h symmetry in rutile, but C 2v in anatase and C 1 in brookite. This implies that U should be different because the d orbitals interact differently. However, to our knowledge, there are no calculated values on U and J for Ti 4+ in TiO 2 . In reproducing the experimental band gap, some authors have found that U = 7.0 eV for Ti 4+ in SrTiO 3 . 91 Other authors report U values but for cases in which Ti has been reduced (by oxygen vacancies 92,93 with GGA/LDA U eff = (U − J) ≥ 4.2 eV or with dopants 94 with LDA U eff = 5.6 eV). It is clear that the use of DFT+U is empirical, but the development of a non-empirical method for assigning the parameter would be useful for taking this forward. In addition, experiments with low temperature calorimetry of the titania phases would be also be valuable for comparison purposes.
In our calculations, the temperature is not taken into account, which can introduce more differences. The calculated structures with GGA/PW91 functional are in agreement with experiments and with previous theoretical work by Beltrán et al. 66 (using B3LYP functional), by Gong and Selloni 87 (using PWSCF functional) by Zhu et al. 95 (using GW methods) and by Mohamad et al. 96 (using Perdew-Burke-Ernzerhof (PBE) method). GGA functionals are known to systematically underestimate the band gap. However, the relative order of the calculated band gaps is also consistent with experiments. We must be cautious about over interpreting this result as it has been shown, for example, that the optical band gap is dependent on particle size and hence the bulk crystal value can differ (anatase: 3.2 eV, 76 3.34 eV). 13 That has led us to the determination of the band gap of surfaces in order to better understand some trends.
The Ti−O bulk distances will play a key role in determining the energetic cost of forming surfaces. Generally, it is energetically more favorable to cut the longest distance, which is given by the surface energy before the ionic relaxation. The (001) surface is the most stable of anatase before relaxation 26 and the (210) surface is the most stable in the case brookite. The relaxation can change the stability order. The relaxation depends on the surface symmetry and the type of Lewis sites. For example, the (001) surface of anatase is flat and has equivalent surface Ti V atoms and equivalent surface µ 2 O atoms. The surfaces with low coordination of Ti species (i.e., Ti IV ) have the higher surface energies. For example, the (301) and (111) surfaces of brookite, having Ti IV sites, are less stable than the (210) surface having only Ti V surface sites.
The surface energies and the coordination of the surface species are therefore linked to their stability and reactivity; less stable is the surface, more reactive it is expected. The comparative adsorption energies of water and formic acid on (101) surface of anatase and (210) surface of brookite confirms this trend. 97 With the surface (210) of brookite being more reactive, it is expected to be a better catalyst compared to an anatase phase exposing (101)  The bulk of the brookite phase has six different distances between 1.873 and 2.081 Å. The variation of the surface Ti−O distances can be used in the interpretation of EXAFS spectra. The main difference between the (301) and (111) surfaces is that the mean Ti−O distances is shorter on the (301) surface and some bonds stretch (ca. 0.1 Å) on the (111) surface which means that the Ti V is predicted to become more acidic. There are also more Ti IV ions on the (301) surface and therefore a higher intrinsic acidity. The brookite offers a wide variety of Ti and O geometries that will result in the different reactivity of each species (more or less acidic and more or less accessible).
As noted above, one visual way of comparing the relative stability of surfaces is to combine the Gibbs-Curie-Wulff law with the calculated surface energies and cell parameters to give the equilibrium morphology of the crystal. In natural samples of rutile, the (110) surfaces are predominant while the common form of anatase is a bipyramid with a squared base. To our knowledge, the brookite phase does not have a typical form but some particles seem to have a large proportion of the (010) surface. 98,99 Ziólkowski 100 did report some typical morphologies of the three phases. The natural crystal form is likely to be more relevant to compare to our calculations because they are more likely grown near thermodynamic conditions. In contrast, the experimentally grown morphology is presumably to be kinetically controlled and hence varies with the synthesis conditions. Thus the nanoparticles may not be at thermodynamic equilibrium. Most of the surfaces have similar surface energies, which might explain why various morphologies have been observed; the morphology depends strongly on the external conditions. However, the prediction of the morphology of the particle at the thermodynamic equilibrium under vacuum gives some insight in order to understand some experimental morphologies. In particular, atomistic simulations enable to calculate a wide range of surfaces and some surfaces not calculated with DFT appear clearly in the morphology. The {101} facet is one of the most stable surface, appearing in the morphology (see Figure 7b) and experimentally observed 21 while the exposed {110} facet (see Figure 7b) having the same surface energy order is observed by Ohno et al. 23 For a direct comparison with different polymorphs, it is crucial to use the same method (for example, on rutile (110), surface energies calculated with different methods vary from 0.3 up to 1.8 J m -2 ). 29,31,32,34,37,51 It is interesting to note that the (110) surface of rutile is slightly less stable than the (101) surface of anatase. This result is surprisingly the opposite of the trend obtained by Lazzeri et al. 29 (anatase (101): 0.44 J m -2 , rutile (110): 0.31 J m -2 ) who have used a PBE functional, whereas their results on anatase surfaces are in complete agreement with a previous study with PW91. 26 The surface energy of the rutile (110) oscillates with the number of layers as it has been demonstrated by other studies. 40 We have calculated the surface energy for different numbers of layers with a large supercell with a surface area of 77.88 Å 2 , k-point meshes: 4 × 3 × 1. The number of Lewis sites is: 8.5 µmol m -2 µ 2 O and 8.5 µmol m -2 Ti V . The number of broken bonds is 16.0 µmol m -2 . For 5 and 6 layers, the surface energy is, respectively, 0.57 and 0.44 J m -2 . On increasing the number of layers, we know that the limit should be between those two values. In reducing the surface area at 19.48 Å 2 in order to be able to increase the number of layers, the surface energy is similar for 5 and 6 layers (respectively, 0. This slab is good enough to have a good description of the energetics and allows us to study the adsorption of small molecules. 58,60 Effectively, for 4 layers, the surface energy is 0.44 J m -2 , so the surface with 3 layers is well converged. This is supported by the experimental observation that the anatase phase is more stable below 14 nm because of the higher surface energy of rutile than anatase. 101,102 The order of the number of Lewis sites might help to rationalize the crystal growth on brookite under different pH and PZC tendency. The order is as following (B = brookite, A = anatase, R = tutile): Order This order is then the same for the number of broken bonds. Following the idea from Hu et al., 80 the most acidic (with lower PZC) should be the anatase (101) surface but with the same order of acidity strength than the rutile (110) surface. If we consider that the anatase phase has an increase of the (001) surface area, in that case, the anatase should have a higher PZC. However, the evaluation of the PZC through a simple relationship with the number of broken bonds does not take into account the hydroxylation state of the surfaces which is the basis of the MUSIC (Multisite Complexation) model. 103 The determination of the hydroxylation can help the interpretation of the PZC of oxides as a function of the shape as described by Jolivet et al. 104 The authors have found a PZC of 5.3 for the anatase but the range of PZC from the literature 105 is wide and is highly dependent on the methodology (titration, electrokinetics), the temperature, the shape of crystal. Globally, it is reported that the PZC values of anatase and rutile are similar which should be true if rutile and anatase crystals are respectively dominated by the (110) and (101) surfaces.
The brookite nanoparticles from Pottier et al. 19 have been synthesized at low temperature using a TiCl 4 precursor with HCl and low pH (pH = 1). The morphology is a lozenge tablet. To obtain this morphology using the Gibbs-Curie-Wulff construction, we suggest that the basal surface is the (001) plane. The lateral surfaces perpendiculars to basal surface are the (210) surfaces ( Figure 9 confirming the experimental morphology from Pottier et al.) 19 In addition, it is possible to construct a morphology respecting Pottier et al. 19 and Ohno et al. 23 morphologies with the following relative surface energies: (001) 1.2, (210) 0.72, (111) 1.08 and (201) 1.08.
Another interesting property derived from electronic structure calculations is the smallest band gap of the different surfaces. With the decrease of the band gap, it becomes easier to form an exciton which is required in dye sensitized solar cell devices. Brookite nanoparticles synthesized by Pottier et al. 19 seem in fact to have the ideal morphology for lithium batteries, exposing mainly the (001) surface giving the best access for lithium diffusion. Indeed, this surface gives a direct access to the open channels for Li intercalation and diffusion along the [001] direction. Similar observations with force field methods have been made on the cathode material, LiFePO 4 where the authors demonstrate that flat particles exposing the (010) surfaces favor the reversible intercalation of Li and diffusion along the [010] open channels. 106 While regarding syntheses of anatase exposing mainly (001) surfaces, such morphologies should be better candidates for photovoltaic systems, exhibiting the most conducting surfaces.
Also, the brookite phase has been found to be a potential photocatalyst, 7,78 which may be explained by the high Lewis acidity of the (010) surface (e.g., 4-fold Ti atoms). However, the OH groups can be unfavorable for catalysis reactions, such as photocatalytic oxidation of trichloroethylene, 107 whereas they will be ideal for photomineralization of toluene. 108

Conclusions
This study describes various surfaces of brookite-TiO 2 which have been experimentally observed where the acid synthesis conditions seem to represent an excellent approach for generating highly reactive surfaces. DFT simulations indicate that the brookite phase is a potential catalyst because the main exposed surfaces have acid sites Ti IV (i.e., (111) surface). The (010) surface of brookite exhibits a high concentration of Ti IV and therefore is expected to be a reactive surface. This surface participates at the good photochemical properties as photoinduced hydrophilicity. 20 In contrast, the main traditional surfaces of rutile and anatase have Ti V surface atoms. However, the number of Lewis sites is not a unique index to describe the acido-basic properties of oxides. The quality of the site also plays a key role. The brookite is more hydrophilic and reactive due to exposed surfaces with Ti IV sites.
The understanding of the crystal growth is crucial to be able to select the best surface for a chosen purpose. Regarding our calculations, in comparing morphologies of brookite nanoparticles, those exhibiting a higher surface area of (001) surfaces are predicted to be more efficient for the anode in lithium batteries. Further experimental investigations should be able to confirm our predictions. In addition, another significant result is that the (001) surface of the anatase phase has the lowest calculated band gap and then we suggest that this surface should be more effective for photocatalysis.

Supplementary Information
Supplementary information on surfaces energies using force fields and DFT+U calculations is available free of charge at http://jbcs.sbq.org.br as PDF file.