ARTICLES This is the peer reviewed version of the following ar3cle: Nature Materials 2017, 16, 328-334, which has been published in final form at hFps://www.nature.com/ar3cles/nmat4804. This ar3cle may be used for non-commercial purposes in accordance with Macmillan Terms PUBLISHED ONLINE: 21 and Condi3ons for Self-Archiving. NOVEMBER 2016 | DOI: 10.1038/NMAT4804 Entropic contributions enhance polarity compensation for CeO2(100) surfaces Marçal Capdevila-Cortada* and Núria López* Surface structure controls the physical and chemical response of materials. Surface polar terminations are appealing because of their unusual properties but they are intrinsically unstable. Several mechanisms, namely metallization, adsorption, and ordered reconstructions, can remove thermodynamic penalties rendering polar surfaces partially stable. Here, for CeO2 (100), we report a complementary stabilization mechanism based on surface disorder that has been unravelled through theoretical simulations that: account for surface energies and configurational entropies; show the importance of the ion distribution degeneracy; and identify low di￿usion barriers between conformations that ensure equilibration. Disordered configurations in oxides might also be further stabilized by preferential adsorption of water. The entropic stabilization term will appear for surfaces with a high number of empty sites, typically achieved when removing part of the ions in a polar termination to make the layer charge zero. Assessing the impact of surface disorder when establishing new structure–activity relationships remains a challenge. S urface terminations are responsible for the chemical and physical functionalities of materials1 . Structure–activity relationships, crucial to materials design2,3 , rely on the detailed characterization of the active surface stoichiometry and local geometry for all the exposed facets of the material. In the last decade, shape-directed synthetic protocols4–6 have enabled the production of tailored nanoparticles with facet control beyond the minimum energy Wul morphology (thermodynamic equilibrium) limitations. These nanostructures7–9 can exhibit high-energy surfaces, many of which may be polar and, in principle, unstable. Polar surfaces present macroscopic polarization along the surface normal, which requires compensation to render them stable10 . The compensation strategies11 are of chemical and physical origin. Chemical procedures include the adsorption of background gases12 , mostly leading to hydroxylated surfaces13 , metal-deposited systems14 , which were identified for ZnO(0001) or MgO(111), and the intrinsic modification of the surface stoichiometry. Physical compensation mechanisms encompass an intrinsic modification of the electronic structure. This is achieved by modifying the covalency in the surface layers, inducing their metallization, and was suggested for Fe3 O4 (100) (ref. 15). Alternatively, this material redistributes subsurface cations to interstitial sites forming subsurface vacancies16 . For polar surface terminations, when the cleavage plane is rich in one of the components, for instance oxygenterminated surfaces, polarity compensation can be retrieved by partially eliminating the excess of this component, re-establishing the stoichiometry. This mechanism is typically accompanied by long-range geometric rearrangements that appear as ordered reconstructions. This results in a variety of ordered structures such as triangular nanoislands and honeycomb patterns for ZnO(0001) (refs 13,17), the octopole reconstruction for MgO(111) (ref. 18), or coexisting (nx1) patterns exposing TiO4 pyramids in SrTiO3 (110) (ref. 19). The last perovskite oxide has been the object of many studies as it presents a myriad of surface reconstructions under various experimental conditions, even on the simple low-index (001) surface20 , where also disordered glass-like structures were observed21 , in resemblance to Potts-type models22 . Interestingly, for its polar (111) surface, configurational disorder was shown to play an important role at very high temperatures (above 1,280 C; ref. 23). Cerium oxide has a fluorite structure that shows charge separation. By cutting the lattice along {100} directions, polar unreconstructed O-terminated (O-t) and Ce-terminated (Ce-t) planes appear, Fig. 1. Originally, the CeO2 (100) oriented films and cleaved single crystals were reported to expose oxygen atoms24–26 . This reconstruction presented a half oxygen monolayer termination (O-t), matching a (2 ⇥ 2)R45 or c(2 ⇥ 2) pattern (Fig. 1). Post-annealing treatment of the films under various atmosphere compositions produced di erent stable oxygen-terminated reconstructions27 . Shape control can drive the formation of CeO2 nanocubes28 exposing the polar (100) facets. Such nanostructures present a unique catalytic activity29–31 . Recent scanning tunnelling microscopy (STM) and high-resolution transmission electron microscopy (HRTEM) imaging of ceria nanocubes has revealed simultaneous regions with either cerium or oxygen terminations32,33 . Given the high energy of the half monolayer Ce-t reconstruction34 , an alternative (2 ⇥ 2) termination composed by capping CeO4 pyramids (grown on the O-t) was suggested (Fig. 1), which is slightly more stable than the half monolayer O-t oxygen termination33 . The experiments thus identify various reconstructions that are highly dependent on the sample history (that is, synthesis, support, annealing) and do not provide conclusive stabilization mechanisms. Here, we provide unprecedented evidence and quantification of the stabilization of polar surfaces derived from surface disorder, which appears as a consequence of the elimination of half of the surface oxygen atoms that ensures charge balance. We will demonstrate that the configurational entropy terms derived from the existence of fewer atoms than available positions are indeed crucial to determine the most stable termination at finite temperatures. To this end, we used state-of-the-art density functional theory to obtain the surface energy of the di erent terminations, including vibrational Institute of Chemical Research of Catalonia (ICIQ), The Barcelona Institute of Science and Technology, Av. Països Catalans 16, 43007 Tarragona, Spain. *e-mail: mcapdevila@iciq.es; nlopez@iciq.es 328 NATURE MATERIALS | VOL 16 | MARCH 2017 | www.nature.com/naturematerials © 2017 Macmillan Publishers Limited, part of Springer Nature. All rights reserved. ARTICLES NATURE MATERIALS DOI: 10.1038/NMAT4804 O-t polar Ce-t polar Reconstruction O-t (2 × 2) CeO4-t (2 × 2) Ce-t (2 × 2) 2.00 Ce-t γ (J m−2) 1.80 1.60 O-t CeO4-t (2 × 2) 1.40 0 5 10 15 Percentage of CeO4 20 25 Figure 1 | The di￿erent reconstructions of CeO2 (100). (Top) Fluorite CeO2 lattice with the two polar unreconstructed O- and Ce- (100) surface cleavages of ceria on the sides (100% O or Ce content). Oxygen atoms are represented as red spheres and Ce as yellow ones. (Centre) From these surface cuts, the previously identified O-t (50% O), Ce-t (50% Ce) and pyramidal CeO4 -t (2 ⇥ 2) (12.1%) reconstructions are shown. (Bottom) Surface energy of the O-t, Ce-t, and of CeO4 -t (nxm) reconstructions as a function of the percentage of CeO4 surface area. Open symbols correspond to plain PBE + U values, whereas filled symbols refer to 0 , from HSE06 calculations including thickness extrapolation and the ZPVE. and configurational entropy contributions. In addition, the easy oxygen di usion observed in the Born–Oppenheimer molecular dynamics (BOMD) simulations ensures that there are no kinetic limitations to the observed surface disorder. To analyse the stability of the di erent facets we computed the surface energies, 0 , of the O-t, Ce-t, and CeO4 -t (2 ⇥ 2) reconstructions at zero kelvin. These values include: the reference HSE06 screened hybrid functional; a thickness extrapolation to minimize the contribution from finite size slabs; and the zeropoint vibrational energy (ZPVE) contribution of the atoms on the surface. Details on the thickness extrapolation are provided in the Supplementary Information. The obtained values are 1.709, 1.916 and 1.634 J m 2 , for O-t, Ce-t and CeO4 -t (2 ⇥ 2) respectively. It is worth noting that the slabs are stoichiometric and therefore the 0 values are independent of the chemical potential or partial pressure of its constituent elements, as already demonstrated in ref. 33. The CeO4 -t reconstruction is based on a (2 ⇥ 2) supercell, as suggested by STM experiments33 . However, as this reconstruction is composed of CeO4 pyramids added to the O-t reconstruction, it can be seen as a hybrid of the O-t and Ce-t reconstructions, where these two terminations are indeed the lower (0%) and upper (24%) limits of CeO4 pyramids. The surface energies of other pyramid densities were also obtained, all of which exhibit very similar surface energy values (Fig. 1). Indeed, the CeO4 -t surface exposing 6% of pyramids (that is, a (4 ⇥ 4) unit cell with two alternated pyramids) is 0.024 J m 2 more stable than the (2 ⇥ 2) CeO4 -t. As the energy di erence between the various CeO4 -t distributions is small, and given that the (2 ⇥ 2) was assigned to the experimental STM observations, the CeO4 -t study will be restricted to the (2 ⇥ 2) cell. The CeO4 -t is therefore the most stable reconstruction at 0 K, only slightly more stable than the O-t, as already reported33 . These values lack vibrational contributions, other than the ZPVE. The vibrational entropy contribution, TSvib , for both reconstructions at 300 K is 0.221 J m 2 for the O-t and 0.217 J m 2 for the CeO4 -t (see Supplementary Information for details). The former termination exhibits slightly higher Svib as its surface oxygens are less coordinated. Overall, at 300 K, the values are 1.487 and 1.417 J m 2 for O-t and CeO4 -t, respectively. The O-t surface possesses soft modes involving surface oxygens, in the range of 140–160 cm 1 . This suggests that they may present NATURE MATERIALS | VOL 16 | MARCH 2017 | www.nature.com/naturematerials © 2017 Macmillan Publishers Limited, part of Springer Nature. All rights reserved. 329 ARTICLES NATURE MATERIALS DOI: 10.1038/NMAT4804 O-t 600 K CeO4-t 800 K 800 K Figure 2 | Surface dynamics of the O-t and CeO4 -t reconstructions. (Top) Two-dimensional probability distribution functions (2D-PDF) of lattice oxygens for the O-t and CeO4 -t reconstructions, obtained from the BOMD simulation at each temperature. Red regions refer to the highest recurrence whereas dark blue translates to no recurrence. The top view of both initial reconstructions is also shown in the background. The 2D-PDF of the CeO4 -t does not include the topmost Ce atom. (Centre) Trajectories of each surface oxygen from the BOMD simulation. (Bottom) Peak amplitudes of each surface oxygen PDF. larger thermal fluctuations than those in CeO4 -t, due to the inherent lower coordination as a result of the surface reconstruction. Indeed, the BOMD simulations on the two reconstructions, at 400, 600, 700 and 800 K show that surface oxygens (Os ) exhibit larger thermal ellipsoids on O-t compared with CeO4 -t (Fig. 2 and Supplementary Fig. 3). Raising the temperature induces an enhancement on the thermal fluctuations in O-t, whereas there is almost no increase in the CeO4 -t. Indeed, the Os of the latter reconstruction exhibit very similar thermal behaviour to the subsurface oxygens of the O-t reconstruction (Supplementary Fig. 4). More importantly, Os di usion appears at 600 K on the O-t reconstruction, and becomes very significant at 800 K (Fig. 2). Similar oxygen and cerium di usion processes were observed using HRTEM on edges and corners of {100} facets on ceria nanoparticles35,36 . On the contrary, the CeO4 -t does not present any Os di usion even at 800 K. These results show that, at commonly explored temperatures, ceria(100) presents di usion of its surface oxygens, and that this process occurs on the well-defined (100) planes being not restricted to edges or singular motifs. The di usion channels of surface oxygens were further investigated by retrieving the thermodynamics, 1E, and kinetics, Ea , for some selected processes. Starting from the O-t structure, the 330 di usion of one Os to an adjacent vacancy site is slightly endothermic, 0.14 eV, and exhibits an energy barrier of 0.45 eV. This Ea value is moderately lower than the energy barrier for the di usion of an oxygen atom on the reduced ceria bulk, 0.60 eV (ref. 37). The di usion of a second Os was also investigated, for the most relevant Os (Supplementary Fig. 6). The di usion barrier values for these subsequent steps present similar or even lower barriers, ranging from 0.23 to 0.87 eV (Supplementary Fig. 6). In turn, the di usion of the two types of surface oxygen on the CeO4 -t reconstruction is more di cult. For the uncoordinated ones (Os,un ), Ea is 0.80 eV on a 0.50 eV endothermic process, whereas Ce-coordinated oxygens Os,Ce require 1.88 eV on a 1.53 eV endothermic process (Supplementary Fig. 6). Therefore, Os di usion on this reconstruction is a more energy-demanding process compared with the O-termination, which explains why di usion is not observed in the BOMD simulations. The configurations derived from the di usions on the O-t surface present very small energy di erences, indicating that some configurations may coexist. Furthermore, these O-t distributions are actually interchangeable, given the low oxygen di usion barriers and as shown by the BOMD simulations. On the basis of Boltzmann statistics, each configuration has a probability of occurrence, Pm , NATURE MATERIALS | VOL 16 | MARCH 2017 | www.nature.com/naturematerials © 2017 Macmillan Publishers Limited, part of Springer Nature. All rights reserved. ARTICLES NATURE MATERIALS DOI: 10.1038/NMAT4804 1.0 1.9 0.8 1.7 1.6 Probability, Pm γ (J m−2) 1.8 0.6 0.4 1.5 1.48 0.2 Distributions Pm < 0.03 1.47 1.46 1.45 0.0 0 200 400 T (K) 600 800 Figure 3 | Surface energy and probability of all O-t (3 ⇥ 3) oxygen distributions. (Left) Surface energies of the 398 non-degenerate (3 ⇥ 3) O-t configurations at the PBE + U level and the resulting density of states once balanced by their degeneracy. The four main configurations are highlighted in red, green, orange and blue. The bottommost panel is a zoom of the lower part of the main panel, shown for clarity. (Right) Probability of occurrence (Pm ) as a function of the temperature for the four main configurations (depicted in the rightmost panel). which is a function of the temperature and depends on the degeneracy of the configuration: ✓ ◆ ⌦m Em Pm = exp (1) Z kB T where Z is the canonical partition function, ⌦m is the degeneracy, and Em is the energy of the configuration m. Since the (4 ⇥ 4) cell has more than 109 configurations, its associated distributions cannot be evaluated quantitatively. The (3 ⇥ 3), on the other hand, presents 48,620 configurations. We used the Site Occupation Disorder algorithm38 to identify the irreducible configurations, that is, not related by symmetry operations, on the O-t (3 ⇥ 3) surface. The approximately 50,000 initial configurations are then reduced to 398 that are unique. The multiplicity of these configurations ranges from 2 to 144. These unique configurations were calculated to take into account their contribution to the surface energy with respect to the temperature. The surface energies of the 398 configurations of the O-t are shown in Fig. 3. To reduce the computational burden, the calculations were performed at the PBE + U level, on thick slabs of 13 surface layers to diminish size e ects. Although PBE + U tends to underestimate surface energy values, it successfully captures the right trend as shown in Fig. 1. The energy distribution for these configurations exhibits almost a continuum spectrum. The lowest energy configuration starts at 1.45 J m 2 , and then there is a large group of almost continuous configurations between 1.52 and 1.60 J m 2 ; this fraction accounts for 148 of the 398 unique configurations. To reproduce the real energy spectrum the contributions are weighted by their multiplicities (Fig. 3). On average, the spectrum features a distribution similar to a Maxwell– Boltzmann distribution, with a maximum at 1.58 J m 2 . The most stable configuration is the regular O-t (in red in Fig. 3), very close in energy to an oxygen zigzag configuration (in green in Fig. 3). The structures following in the stability order also present similar zigzag patterns. The probabilities of the most stable configurations with respect to the temperature are shown in Fig. 3. The regular O-t is the most stable reconstruction at 0 K but its degeneracy is only 2. The aforementioned zigzag configurations present higher multiplicities and therefore become more likely at higher temperatures; in fact, at 600 K the probability of the regular O-t is only 0.05. The second most stable configuration has ⌦ = 72 and hence it becomes the most probable distribution at 190 K. Its probability reaches a maximum at 395 K (0.87) but many other configurations with the same or higher degeneracy achieve probabilities between 0.01 and 0.05 at high temperature and decrease that of the main configuration. Then, the surface energy obtained for the O-t surface needs to be re-evaluated to account for the contributions coming from the configurations that might be simultaneously present for this stoichiometry. These additional stabilizing contributions stem from the configurational entropy (Sconf ) terms present on the O-t surface. The configurational entropy related to each distribution is Sconf,m = kB ln(⌦m ), and the total configurational entropy becomes Sconf = 6m Pm Sconf,m (further details are given in the Methods and the Supplementary Information). The surface energy can be now expressed as a function of the temperature (Fig. 4), by means of the equation (T ) = 0 (T /A)(Svib + Sconf ). The entropic contributions then stabilize the O-t over the CeO4 -t reconstruction, and even though the latter is more stable at 0 K, their relative stability decreases when the temperature increases until the O-t becomes more stable at 790 K. Considering the intrinsic errors in the current implementations of density functional theory39 , we have identified that two structures present the same surface energy when their surface energy di erence is lower than |1 | < 0.001 eV Å 2 . Under this premise, both the O-t and CeO4 -t reconstructions are isoenergetic on the 670–905 K temperature range (gradient region of Fig. 4a), and the O-t encompasses various oxygen configurations that are interchangeable. In comparison, oxygen loss occurs at 963 K for pO2 = 1 bar. The di erent chemical response of the O-t and CeO4 -t surfaces might also have an important e ect when examining the surface termination and can explain the dependence on the history of the sample in the observed surface termination. Water is the clearest example of a probe molecule, as it is typically present as a residual impurity in the ultrahigh-vacuum chambers where the most detailed STM experiments are carried out. Thus, certain water content in the samples is unavoidable, even under ultrahigh-vacuum conditions40 . Water adsorbs dissociatively on both reconstructions, with higher adsorption energy on the regular O-t ( 1.60 eV) than on the CeO4 -t ( 1.39 eV). This fact translates into a higher stabilization of the O-t reconstruction compared NATURE MATERIALS | VOL 16 | MARCH 2017 | www.nature.com/naturematerials © 2017 Macmillan Publishers Limited, part of Springer Nature. All rights reserved. 331 ARTICLES a NATURE MATERIALS DOI: 10.1038/NMAT4804 2.00 γ (J m−2) 1.50 1.00 0.50 0.00 O-t CeO4-t 0 200 400 600 800 1,000 T (K) O-t −12 −20 0.0 0 c 500 T (K) CeO4-t −4 1.0 θ ln(p/p0) −4 ln(p/p0) b −12 −20 1,000 0 500 T (K) 0.04 −4 O-t O-t 0.02 −8 0.00 −12 −0.02 −16 −0.04 CeO4-t −20 γ CeO4-t − γ O-t (J m−2) ln(p/p0) 1,000 0 200 CeO4-t 400 600 800 1,000 −0.06 T (K) Figure 4 | Surface energy including entropy contributions and water adsorption. a, Surface energies for the O-t (light blue) and CeO4 -t (brown) reconstructions as a function of the temperature. The ZPVE and entropic (Svib and Sconf ) terms are included. The gradient area shows the temperature range where both reconstructions are almost isoenergetic 2 (|1 | < 0.001 eV Å ). b, Contour plot of the water coverage (✓ ) as a function of the water pressure, ln (p/p0 ), and the temperature, T. c, Contour O-t ) as a function of the plot of the surface energy di￿erence ( CeO4 -t water pressure, ln (p/p0 ), and the temperature, T. Negative values (red–purple) indicate more stable CeO4 -t reconstruction, whereas positive values (orange–green) point to the region where O-t is thermodynamically more stable. with the CeO4 -t following water interaction. Therefore, an extra stabilizing contribution can be introduced in the surface energy, the ads term, defined as the adsorption energy of water per surface area. Thus, the total surface energy can be expressed as T , pH2 O = surf (T ) + ✓ T , pH2 O ads (T , pH2 O ), where surf is the surface energy when no water is present and ✓ is the water coverage41 . Based on a Langmuir adsorption model, the water coverage depends on the water gas pressure, the temperature, and its adsorption energy (further details are given in the Supplementary Information). Since we are interested in simulating low water presence in the samples, the range of pressures was restricted to 332 1.8 ⇥ 10 2 and 2.1 ⇥ 10 9 bar. Within this simple model, water coverage is 1 at low temperatures, for both reconstructions (Fig. 4b). Due to its stronger bond to the O-t, the water coverage remains high for a larger window of temperatures and pressures on this termination. The surface energy of both reconstructions can then be expressed by means of ab initio thermodynamics as a function of the water chemical potential (see Supplementary Information). The di erence in surface energy of the CeO4 -t and O-t reconstructions, CeO4 -t – O-t , is plotted in Fig. 4c as a function of the water pressure, ln(p/p0 ), and the temperature, T . The complex behaviour found is a combination of two factors: the surface energies following water adsorption ads (Supplementary Fig. 7) and the water coverages ✓ (Fig. 4b). The ads term is more stable for the O-t than the CeO4 -t, for all T and p ranges, while the coverage ✓ is higher for O-t in a wider window than for CeO4 -t. Thus, for intermediate pressures (ln(p/p0 ) = 8.0) and low temperatures, where the water coverage is high on both CeO4 -t and O-t, surf is the dominant term for CeO4 -t to become more stable. Starting at 350 K, the di erence in ads becomes larger than the di erence in surf , thus prompting the O-t to become more stable. In addition, the water coverage in CeO4 -t turns zero at a lower temperature than in O-t, which increases the stability of the O-t in the range where ✓CeO4 -t = 0 and ✓O-t is still high. Then, when water coverage is zero for both reconstructions, at around 600 K, CeO4 -t becomes more stable again, until at about 790 K, where the entropy terms invert the stability order again and we retrieve the crossing found in Fig. 4a. Severe discrepancies have been found in the microscopy experiments of the surface of CeO2 (100) due to the di culties in imaging and the di erent sample history. It was suggested that ceria(100) exhibits several terminations because of their small energy di erences and that they may be in constant flux42 . Now we have proved that these di erences have a physical origin (attributable to configurational entropies) that can be further enhanced by the di erential chemical properties of the two surfaces. To facilitate further characterization of these surfaces belonging to O-t configurations, STM and HRTEM simulations of those with probability higher than 0.05 at 500 K are provided in Fig. 5. The aforementioned zigzag patterns are clearly identified in the STM images and they were also found after di usion at 800 K in the BOMD simulations on the (4 ⇥ 4) supercells. Any of these compensation mechanisms modifies the surface and thus perturbs the chemical and physical responses. The ultimate consequence of the pseudorandom distribution of the anions on the O-t reconstruction is that the local configuration of the reactive sites might control the adsorption modes for multiple functionalized molecules, changing the adsorption configuration, and reactive gases can modify the surface local configurations as the energies involved are very low. The physical properties will also be similarly a ected, as the CeO4 -t shows lower unoccupied states that belong to the exposed Ce (see Supplementary Fig. 7). Thus, the chemistry of ceria nanocubes can be a ected by all these contributions. The important contribution of the surface configurational entropy in the surface termination can be applied to other materials, particularly oxides. There are two main thermodynamic requirements for this behaviour: first, that the anion-polar termination can be stabilized by the loss of some terminal anions; as a consequence, the number of positions is larger than the number of particles, which ensures a favourable Sconf ; and, second, the energy di erence between these configurations must be small. Since the sites are almost equivalent, in principle this requirement should be possible; however, dense areas with high anion concentration (or, analogously, with high depletion) might be high in energy. There is a further kinetic constraint that might limit the equilibration between the di erent configurations; as a consequence, the di usion of the particles shall allow the mass transport required for the configurations to be present. These results can be generalized to NATURE MATERIALS | VOL 16 | MARCH 2017 | www.nature.com/naturematerials © 2017 Macmillan Publishers Limited, part of Springer Nature. All rights reserved. ARTICLES NATURE MATERIALS DOI: 10.1038/NMAT4804 x y z x Figure 5 | Surface characterization of the most stable O-t distributions. Simulated STM (top) and HRTEM (bottom) images of the O-t configurations with higher probability at 500 K. other polar surfaces that contain a large number of defects with an even larger number of potentially equivalent sites on open surfaces. Similarly, the behaviour can appear in the context of surface oxide formation. For instance, copper (100) oxidizes forming chains of di erent lengths with local order in the formation of the Cu–O–Cu–O patterns43 but these chains are disordered on the surface. Again, the configurational entropy stems from the number of equivalent sites, which are much larger than the number of Cu–O units. This might make the structure more stable than a sequential oxygen incorporation with the formation of ordered patterns. For compounds such as SrTiO3 , disordered reconstructions were also identified but at a much higher temperature (1,280 C) (ref. 23). Finally, it is worth stressing that although entropy contributions have often been neglected, various studies have recently shown their pivotal role in complex materials44,45 . The new findings are closely connected with the pioneering Potts models for the description of spin glasses or adsorption processes22,46 . In summary, our results indicate that the way open polar surfaces are understood needs to account for the intrinsic disorder. This configurational contribution arises from almost degenerate formation enthalpies of a wide variety of configurations due to the partial removal of one of the components, particularly when the configurations can be equilibrated through fast di usion kinetic processes. The disordered configurations will have a substantial impact on how these materials behave and work when used in electronics or as catalysts. Methods Methods, including statements of data availability and any associated accession codes and references, are available in the online version of this paper. Received 24 April 2016; accepted 18 October 2016; published online 21 November 2016 References 1. Wu, Z. & Overbury, S. H. Catalysis by Materials with Well-Defined Structures (Elsevier, 2015). 2. Nørskov, J. K., Bligaard, T., Rossmeisl, J. & Christensen, C. H. Towards the computational design of solid catalysts. Nat. Chem. 1, 37–46 (2009). 3. Calle-Vallejo, F., Lo reda, D., Koper, M. T. M. & Sautet, P. Introducing structural sensitivity into adsorption–energy scaling relations by means of coordination numbers. Nat. Chem. 7, 403–410 (2015). 4. Sun, Y. & Xia, Y. Shape-controlled synthesis of gold and silver nanoparticles. Science 298, 2176–2179 (2002). 5. Xia, X. et al. On the role of surface di usion in determining the shape or morphology of noble-metal nanocrystals. Proc. Natl Acad. Sci. USA 110, 6669–6673 (2013). 6. Lu, M. et al. Shape-controlled synthesis of hybrid nanomaterials via three-dimensional hydrodynamic focusing. ACS Nano 8, 10026–10034 (2014). 7. Xia, Y., Xiong, Y., Lim, B. & Skrabalak, S. E. Shape-controlled synthesis of metal nanocrystals: simple chemistry meets complex physics? Angew. Chem. Int. Ed. 48, 60–103 (2009). 8. Wei, Z. & Matsui, H. Rational strategy for shaped nanomaterial synthesis in reverse micelle reactors. Nat. Commun. 5, 3870 (2014). 9. Helmi, S., Ziegler, C., Kauert, D. J. & Seidel, R. Shape-controlled synthesis of gold nanostructures using DNA origami molds. Nano Lett. 14, 6693–6698 (2014). 10. Tasker, P. W. The stability of ionic crystal surfaces. J. Phys. C Solid State Phys. 12, 4977–4984 (1979). 11. Noguera, C. & Goniakowski, J. Polarity in oxide nano-objects. Chem. Rev. 113, 4073–4105 (2013). 12. äodziana, Z., Topsøe, N.-Y. & Nørskov, J. K. A negative surface energy for alumina. Nat. Mater. 3, 289–293 (2004). 13. Lauritsen, J. V. et al. Stabilization principles for polar surfaces of ZnO. ACS Nano 5, 5987–5994 (2011). 14. Goniakowski, J. & Noguera, C. Characteristics of Pd deposition on the MgO(111) surface. Phys. Rev. B 60, 16120–16128 (1999). 15. Fonin, M. et al. Surface electronic structure of the Fe3 O4 (100): evidence of a half-metal to metal transition. Phys. Rev. B 72, 104436 (2005). 16. Bliem, R. et al. Subsurface cation vacancy stabilization of the magnetite (001) surface. Science 346, 1215–1218 (2014). 17. Dulub, O., Diebold, U. & Kresse, G. Novel stabilization mechanism on polar surfaces: ZnO(0001)-Zn. Phys. Rev. Lett. 90, 016102 (2003). 18. Barbier, A. et al. Atomic structure of the polar NiO(111)-p(2 ⇥ 2) surface. Phys. Rev. Lett. 84, 2897–2900 (2000). 19. Enterkin, J. A. et al. A homologous series of structures on the surface of SrTiO3 (110). Nat. Mater. 9, 245–248 (2010). 20. Erdman, N. et al. The structure and chemistry of the TiO2 -rich surface of SrTiO3 (001). Nature 419, 55–58 (2002). 21. Kienzle, D. M., Becerra-Toledo, A. E. & Marks, L. D. Vacant-site octahedral tilings on SrTiO3 (001), the (13 ⇥ 13) R 33.7 surface, and related structures. Phys. Rev. Lett. 106, 176102 (2011). 22. Kirkpatrick, T. R. & Wolynes, P. G. Stable and metastable states in mean-field Potts and structural glasses. Phys. Rev. B 36, 8552–8564 (1987). 23. Marks, L. D., Chiaramonti, A. N., Rahman, S. U. & Castell, M. R. Transition from order to configurational disorder for surface reconstructions on SrTiO3 (111). Phys. Rev. Lett. 114, 226101 (2015). 24. Herman, G. S. Surface structure determination of CeO2 (001) by angle-resolved mass spectroscopy of recoiled ions. Phys. Rev. B 59, 14899–14902 (1999). 25. Kim, Y. J. et al. Growth and structure of epitaxial CeO2 by oxygen-plasma-assisted molecular beam epitaxy. J. Vac. Sci. Technol. 17, 926–935 (1999). 26. Nörenberg, H. & Harding, J. H. The surface structure of CeO2 (001) single crystals studied by elevated temperature STM. Surf. Sci. 477, 17–24 (2001). 27. Solovyov, V. F. et al. Highly e cient solid state catalysis by reconstructed (001) ceria surface. Sci. Rep. 4, 4627 (2014). 28. Qiao, Z.-A., Wu, Z. & Dai, S. Shape-controlled ceria-based nanostructures for catalysis applications. ChemSusChem 6, 1821–1833 (2013). 29. Vilé, G., Colussi, S., Krumeich, F., Trovarelli, A. & Pérez-Ramírez, J. Opposite face sensitivity of CeO2 in hydrogenation and oxidation catalysis. Angew. Chem. Int. Ed. 53, 12069–12072 (2014). 30. Agarwal, S., Mojet, B. L., Le erts, L. & Datye, A. K. Catalysis by Materials with Well-Defined Structures 31–70 (Elsevier, 2015). 31. Mann, A. K. P., Wu, Z. & Overbury, S. H. Catalysis by Materials with Well-Defined Structures 71–97 (Elsevier, 2015). NATURE MATERIALS | VOL 16 | MARCH 2017 | www.nature.com/naturematerials © 2017 Macmillan Publishers Limited, part of Springer Nature. All rights reserved. 333 ARTICLES NATURE MATERIALS DOI: 10.1038/NMAT4804 32. Lin, Y., Wu, Z., Wen, J., Poeppelmeier, K. R. & Marks, L. D. Imaging the atomic surface structures of CeO2 nanoparticles. Nano Lett. 14, 191–196 (2014). 33. Pan, Y. et al. Ceria nanocrystals exposing wide (100) facets: structure and polarity compensation. Adv. Mater. Interfaces 1, 1400404 (2014). 34. Baudin, M., Wójcik, M. & Hermansson, K. Dynamics, structure and energetics of the (111), (011) and (001) surfaces of ceria. Surf. Sci. 468, 51–61 (2000). 35. Möbus, G. et al. Dynamics of polar surfaces on ceria nanoparticles observed in situ with single-atom resolution. Adv. Funct. Mater. 21, 1971–1976 (2011). 36. Bhatta, U. M. et al. Cationic surface reconstructions on cerium oxide nanocrystals: an aberration-corrected HRTEM study. ACS Nano 6, 421–430 (2012). 37. Dholabhai, P. P., Adams, J. B., Crozier, P. & Sharma, R. Oxygen vacancy migration in ceria and Pr-doped ceria: a DFT+U study. J. Chem. Phys. 132, 94104 (2010). 38. Grau-Crespo, R., Hamad, S., Catlow, C. R. A. & de Leeuw, N. H. Symmetry-adapted configurational modelling of fractional site occupancy in solids. J. Phys. Condens. Matter 19, 256201 (2007). 39. Lejaeghere, K. et al. Reproducibility in density functional theory calculations of solids. Science 351, aad3000 (2016). 40. Li, S.-C., Losovyj, Y., Paliwal, V. K. & Diebold, U. Photoemission study of azobenzene and aniline adsorbed on TiO2 anatase (101) and rutile (110) surfaces. J. Phys. Chem. C 115, 10173–10179 (2011). 41. äodziana, Z., Nørskov, J. K. & Stoltze, P. The stability of the hydroxylated (0001) surface of ↵-Al2 O3 . J. Chem. Phys. 118, 11179–11188 (2003). 42. Mullins, D. R. The surface chemistry of cerium oxide. Surf. Sci. Rep. 70, 42–85 (2015). 43. Gattinoni, C. & Michaelides, A. Atomistic details of oxide surfaces and surface oxidation: the example of copper and its oxides. Surf. Sci. Rep. 70, 424–447 (2015). 334 44. Gludovatz, B. et al. A fracture-resistant high-entropy alloy for cryogenic applications. Science 345, 1153–1158 (2014). 45. Butler, K. T., Walsh, A., Cheetham, A. K. & Kieslich, G. Organised chaos: entropy in hybrid inorganic–organic systems and other materials. Chem. Sci. 7, 6316–6324 (2016). 46. Za™uska-Kotur, M. A. The kinetic Potts model in the description of surface dynamics. Surf. Sci. 265, 196–208 (1992). Acknowledgements This research has been supported by the ERC Starting Grant and Proof of Concepts (ERC-2010-StG-258406, ERC-2015-PoC_680900), the Ministerio de Economía y Competitividad—MINECO (CTQ2015-68770-R), and the Generalitat de Catalunya—AGAUR (SGR-2014-SGR-145). M.C.-C. acknowledges MINECO for a ‘Juan de la Cierva—Formación’ fellowship (FJCI-2014-20568). We acknowledge BSC-RES and CSUC for providing generous computational resources. We thank R. Grau-Crespo (University of Reading) for kindly providing us with the SOD software and A. Selloni (Princeton University) for critically reading the manuscript. Author contributions M.C.-C. performed the calculations. M.C.-C. and N.L. analysed the data and prepared the manuscript. Additional information Supplementary information is available in the online version of the paper. Reprints and permissions information is available online at www.nature.com/reprints. Correspondence and requests for materials should be addressed to M.C.-C. or N.L. Competing financial interests The authors declare no competing financial interests. NATURE MATERIALS | VOL 16 | MARCH 2017 | www.nature.com/naturematerials © 2017 Macmillan Publishers Limited, part of Springer Nature. All rights reserved. ARTICLES NATURE MATERIALS DOI: 10.1038/NMAT4804 Methods All the calculations were performed using periodic boundary conditions at the density functional theory level with the Vienna ab initio simulation package (VASP, version 5.3.3; refs 47,48). The Perdew–Burke–Ernzerhof (PBE)49 GGA functional and the Heyd–Scuseria–Ernzerhof (HSE06)50 screened Coulomb hybrid functional were used. For the PBE calculations, the self-interaction error was diminished by the addition of an e ective Hubbard U term, Ue , to the Ce(4f) states, following the approach of ref. 51, setting the value to 4.5 eV (ref. 52). Projector-augmented wave pseudopotentials53 were used to describe the core electrons with a plane-wave cuto energy of 500 eV for the valence electrons. The electronic and geometry optimization convergence criteria were set to 10 6 eV and 0.015 eV Å 1 , respectively. The transition states for oxygen di usion were located by means of the climbing image nudged elastic band method54 . The lattice parameter was optimized using a denser 0-centred 7 ⇥ 7 ⇥ 7 k-point mesh and a plane-wave cuto value of 700 eV that led to a value of acalc = 5.497 Å, in good agreement with the experimental value (aexp = 5.411 Å) (ref. 55). The O-t, CeO4 -t and Ce-t reconstructions of the (100) surface were initially modelled as (2 ⇥ 2) periodically repeated symmetric slabs separated by 15 Å of vacuum space, where one half of the slab was allowed to relax. Additional (3 ⇥ 3) and (4 ⇥ 4) unit cell calculations were performed for the O-t reconstruction. (3 ⇥ 2), (3 ⇥ 3), (4 ⇥ 2), (4 ⇥ 3) and (4 ⇥ 4) unit cells containing one CeO4 pyramid, and the (4 ⇥ 4) containing two alternated pyramids were also evaluated for the CeO4 -t. The Brillouin zone was described using a 0-centred n1 ⇥ n2 ⇥ 1 k-point mesh, where n1 and n2 were set to 3, 2 or 1 for unit cells with axis 2, 3 and 4, respectively, in the corresponding direction. Born–Oppenheimer molecular dynamics (BOMD) were performed on the (4 ⇥ 4) O-t and (2 ⇥ 2) CeO4 -t reconstructions using the PBE + U functional, under the NVT ensemble. The temperature was set to 400, 600, 700 and 800 K, controlled by a Nosé–Hoover thermostat56,57 . All equilibration runs were performed over 1 ps using a time step of 1 fs. The production runs vary with the temperature: 50, 30, 20 and 15 ps, using a time step of 5, 3, 2 and 1 fs, for 400, 600, 700 and 800 K, respectively. The electronic convergence criterion was increased to 10 7 eV for all BOMD simulations. The evaluation of the configurational entropy and the initial set of distributions were performed using the Site Occupancy Disorder (SOD) program38 . Further details on how the surface energy and the entropy contributions were obtained are given in the Supplementary Information. STM simulations were performed under the Terso –Hamann approximation with constant current58 , at the 10 6 e Å 3 isodensity and simulated bias voltage of 1.0 V (that is, evaluating the density of the filled states from 1.0 eV below the Fermi level to the Fermi level, as we want to focus on the oxygen atoms). The isodensity value used translates into an average distance of about 3 Å to the surface, for the four simulated STM images of Fig. 5. HRTEM images were obtained using the clTEM software with an acceleration voltage of 100 kV, Cs = 1 µm, and 30 nm defocus. The main structural models have been uploaded to the ioChem-BD database59 . Data availability. The data that support the plots within this paper and other findings of this study are available from the corresponding authors on request. The main structural models can be obtained from the ioChem-BD database at http://dx.doi.org/10.19061/iochem-bd-1-12. References 47. Kresse, G. E cient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. Phys. Rev. B 54, 11169–11186 (1996). 48. Kresse, G. & Furthmüller, J. E ciency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set. Comput. Mater. Sci. 6, 15–50 (1996). 49. Perdew, J. P., Burke, K. & Ernzerhof, M. Generalized gradient approximation made simple. Phys. Rev. Lett. 77, 3865–3868 (1996). 50. Krukau, A. V., Vydrov, O. A., Izmaylov, A. F. & Scuseria, G. E. Influence of the exchange screening parameter on the performance of screened hybrid functionals. J. Chem. Phys. 125, 224106 (2006). 51. Dudarev, S. L., Savrasov, S. Y., Humphreys, C. J. & Sutton, A. P. Electron-energy-loss spectra and the structural stability of nickel oxide: an LSDA+U study. Phys. Rev. B 57, 1505–1509 (1998). 52. Fabris, S., de Gironcoli, S., Baroni, S., Vicario, G. & Balducci, G. Reply to ‘‘Comment on ‘Taming multiple valency with density functionals: a case study of defective ceria’’’. Phys. Rev. B 72, 237102 (2005). 53. Blöchl, P. E. Projector augmented-wave method. Phys. Rev. B 50, 17953–17979 (1994). 54. Henkelman, G., Uberuaga, B. P. & Jónsson, H. A climbing image nudged elastic band method for finding saddle points and minimum energy paths. J. Chem. Phys. 113, 9901–9904 (2000). 55. Kümmerle, E. & Heger, G. The structures of C–Ce2 O3+ , Ce7 O12 , and Ce11 O20 . J. Solid State Chem. 147, 485–500 (1999). 56. Nosé, S. A molecular dynamics method for simulations in the canonical ensemble. Mol. Phys. 52, 255–268 (1984). 57. Hoover, W. Canonical dynamics: equilibrium phase-space distributions. Phys. Rev. A 31, 1695–1697 (1985). 58. Terso , J. & Hamann, D. R. Theory and application for the scanning tunneling microscope. Phys. Rev. Lett. 50, 1998–2001 (1983). 59. Álvarez-Moreno, M. et al. Managing the computational chemistry big data problem: the ioChem-BD platform. J. Chem. Inf. Model. 55, 95–103 (2015). NATURE MATERIALS | www.nature.com/naturematerials © 2017 Macmillan Publishers Limited, part of Springer Nature. All rights reserved.