Polaron-Driven Surface Reconstructions

Michele Reticcioli, Martin Setvin, Xianfeng Hao, Peter Flauger, Georg Kresse, Michael Schmid, Ulrike Diebold, and Cesare Franchini University of Vienna, Faculty of Physics and Center for Computational Materials Science, A-1090 Vienna, Austria Institute of Applied Physics, Technische Universität Wien, A-1040 Vienna, Austria Key Laboratory of Applied Chemistry, Yanshan University, Qinhuangdao 066004, People’s Republic of China (Received 24 May 2017; revised manuscript received 7 July 2017; published 25 September 2017)

Unraveling the physical mechanism that drives such surface reconstructions is essential to understand physical and chemical properties of surfaces, and it helps optimize materials performance in applications such as microelectronics, catalysis, fuel cells, or surface engineering.The stabilization process typically involves the displacements of surface atoms and electric charge transfer between surface atoms, and often results in nonstoichiometric reconstructions; this mechanism can be affected by defects and heating treatments, and by the interaction of orbitals, lattice, and spins [1,12].In this work, we find that charge trapping can cause surface reconstructions as well.
Excess charge in a flexible and ionic lattice is often present in the form of small polarons, i.e., trapped electrons or holes coupled with the lattice phonon field.Manifestations include a local alteration of the bond lengths, a change of the formal valence at the specific polaronic site, and the emergence of a characteristic peak well localized in the gap region [13][14][15][16][17]. Polarons play a decisive role in electron transport, optical absorption, and chemical reactivity, and have crucial implications in other diverse phenomena including high-T c superconductivity [18], colossal magnetoresistance [19,20], thermoelectricity [21], and photochemistry [22].Formation of polarons is particularly favorable in transition-metal oxides, owing to the strength of the electron-phonon interaction, and is further promoted in the vicinity of a surface, where the crystal lattice is more flexible and the necessary lattice relaxations cost less energy [23][24][25].
Titanium dioxide, TiO 2 , is an archetypal polaron system.The surfaces of this material have been studied extensively in both theory and experiment [24,26,27].TiO 2 single crystals are often reduced, either chemically (by removing O from the lattice, i.e., by creating an oxygen vacancy V O ) or electronically (by adding excess electrons through the introduction of extrinsic dopants).This reduction produces polarons that typically reside in subsurface Ti 3þ sites, although hopping can be thermally activated [17,24].
For highly reduced rutile TiO 2 ð110Þ samples, the surface assumes a (1 × 2) reconstruction rather than a (1 × 1) termination that is observed on more stoichiometric crystals [4,28,29].This (1 × 2) structure is then the low-temperature phase; it can reversibly switch to the (1 × 1) phase when the sample temperature is raised above approximately 1000 K [30], and the transition is facilitated by the surface-to-bulk mass flow [31].Here, we evaluate the surface structure as a function of oxygen vacancy concentration (at the surface), c V O , using density functional theory (DFT), first-principles molecular dynamics (FPMD), scanning tunneling microscopy (STM), and atomic force microscopy (AFM).We find that, with increasing V O density, the resulting polarons tend to form an optimum 3 × 1 arrangement in the (1 × 1) surface.When c V O is further raised, the repulsive interaction between these negatively charged quasiparticles increases the surface free energy.The resulting energy instability is overtaken by the geometrical reconstruction into the (1 × 2) phase.

A. Ab initio calculations
We used DFT and FPMD to investigate the (1 × 1) and (1 × 2) surface phases of the TiO 2 ð110Þ rutile.The calculations were performed using the Vienna ab initio simulation package (VASP) [32,33].We adopted the generalized gradient approximation (GGA) within the Perdew, Burke, and Ernzerhof parametrization [34], with the inclusion of an on-site effective U of 3.9 eV on Ti d-states.In our previous work [17], we calculated this value within the constrained random phase approximation (cRPA) [35] and confirmed that it correctly predicts the experimentally observed degree of localization in the TiO 2 polymorphs rutile and anatase.We also checked the robustness of our conclusions against other values of U (see Ref. [36]).The rutile TiO 2 ð110Þ surfaces were modeled with an asymmetric slab containing five Ti layers with a large twodimensional (2D) 9 × 2 unit cell, resulting in more than 500 atoms.The bottom two layers were kept fixed, whereas all the other atomic sites were relaxed using standard convergence criteria with a plane-wave energy cutoff of 300 eV (lowered to 250 eV for the FPMD runs) and using the Γ point only for the integration in the reciprocal space.For the (1 × 1) phase, up to nine oxygen vacancies (V O 's) were homogeneously included in the top layer (S0) at different concentrations (c V O ¼ 5.6%, 11.1%, 16.7%, 22.2%, 33.3%, 38.9%, 50.0%).This V O distribution was chosen by taking into account the experimentally determined, spatially resolved, autocorrelation functions.More details on this choice are provided in Ref. [36].Each V O supplies two excess electrons, which form two polarons.The (1 × 2) phase was constructed according to the Ti 2 O 3 model [4,37,38] by placing the reconstructed row on top of the five layers, which gives rise to two polarons per (1 × 2) unit cell.
The stability of the surface, and the thermally activated hopping behavior of the polarons, were inspected by combining FPMD at T ¼ 700 K and static T ¼ 0 K DFT calculations using the following strategy: First, for each phase and each c V O , multiple series of FPMD runs were conducted, each starting from an arbitrary polaronic configuration and lasting 10 ps with a time step of 1 fs.This leads to 10 000 different high-temperature structures for each FPMD run.Second, we selected all inequivalent polaronic configurations from this set; about 200 for each phase and V O concentration.Finally, we relaxed these at T ¼ 0 K and found the energetically most stable structure by comparing the total DFT energies.This procedure was repeated for the nine different V O concentrations and for the Ti 2 O 3 reconstructed model, by using the 9 × 2 cell.A slightly lower ground-state energy for the c V O ¼ 16.7% case was obtained by adopting a 6 × 2 unit cell, which allows for a more homogeneous distribution of oxygen vacancies as compared to the 9 × 2 cell at this particular c V O .The charge densities of polarons in the most stable configurations were used for simulating STM images within the Tersoff-Hamann approach [39].The resulting isosurface representations allow for comparison with our experimental data.
The total energies of the most stable configurations were used to build a global surface phase diagram using standard ab initio atomistic thermodynamics [40].By neglecting configurational entropy and phonon contributions, the surface free energy ΔG can be approximated as follows: where E loc is the DFT total energy at T ¼ 0 K for the surface slab with trapped polarons, and E bulk TiO 2 is the energy for the TiO 2 bulk.The adimensional prefactor 1=A scales the energy to the (1 × 1) surface cell.The weighted difference between the number of Ti (n Ti ) and O (n O ) atoms accounts for the deviation from the stoichiometric formula.The chemical potential μ O of oxygen atoms was considered in terms of its deviation Δμ O from the total energy E O 2 of an isolated oxygen molecule: We also constructed an additional phase diagram describing the stability of TiO 2 ð110Þ surfaces forced to have all the electrons delocalized.The delocalized solution was achieved by performing non-spin-polarized calculations.
The surface free energy for the systems with all electrons delocalized is given by replacing the E loc term in Eq. ( 1) with the so-obtained, delocalized, ground-state energy E deloc relax , where the atoms were allowed to relax again.
Moreover, the most stable configurations were analyzed in terms of the polaron formation energy E POL .E POL represents the convenience for the system to trap electrons, and it is calculated as The stability of a polaron, quantified by E POL , is the result of the competition between the structural cost needed to distort the lattice in order to accommodate polarons (E ST ) and the electronic energy gained by localizing the electron in the distorted lattice (E EL ): where E ST is defined as where E deloc constr is the energy of the system forced to have only delocalized electrons and constrained into the structure of the system hosting polarons.
Finally, the results of the FPMD runs were also used to perform a statistical analysis of the polarons and their hopping behavior.For each phase and c V O , we determined the number of occurrences of charge trapping in each layer.We computed the polaron-polaron correlation function, defined as the distribution of the distance between two polarons at a given time step, averaged over all FPMD steps.

B. Experiments
Combined STM/AFM measurements were performed at T ¼ 78 K in an ultrahigh vacuum (UHV) chamber with a base pressure below 2 × 10 −9 Pa, equipped with a commercial Omicron q-Plus LT head.Synthetic single crystal samples (Crystec) were cleaned in an adjoining preparation chamber by repeated sputtering with 1 keV Ar þ ions and annealing to 1000 K.All data were measured on a sample that was strongly reduced by more than 100 cycles of prior sputtering and annealing.Tuning-fork-based AFM sensors with a separate wire for the tunneling current were used (k ¼ 3750 N=m, f 0 ¼ 47 500 Hz, Q ≈ 50 000) [41].Glued to each tuning fork was an electrochemically etched W tip, which was cleaned in situ by field emission and selfsputtering in 10 −4 Pa argon [42].The tips were purposely functionalized by an O adatom.Here, the rutile surface was exposed to O 2 , and bias pulses applied above adsorbed O 2 molecules occasionally resulted in a single O atom being picked up.Comparing force-distance curves (see Ref. [36]) measured above the clean (1 × 1) surface to previous experimental AFM works [43,44] leads us to conclude that the tip is terminated by a single-coordinated O atom (O adatom).In our experience, this is the only tip model where the attractive force above the highest O atoms (O 2c ) is negligible and the repulsive regime is entered directly.With CO-functionalized tips [45], the same AFM contrast was obtained above the (1 × 2) reconstruction, but no good images were obtained on the (1 × 1) phase.

III. RESULTS
The transition from a (1 × 1) to a (1 × 2) structure at the rutile TiO 2 ð110Þ surface is well documented in the literature [4,28].The large-area STM image in Fig. 1 shows a surface typical for a strongly reduced TiO 2 ð110Þ sample.The bulkterminated (1 × 1) structure with isolated oxygen vacancies (V O 's) coexists with stripes of the (1 × 2) reconstruction.When a pristine (stoichiometric) TiO 2 sample is introduced into UHV, it shows only the (1 × 1) surface.The standard sample treatments described above, i.e., sputtering and high-temperature annealing, reduce the crystals and, at first, result in the appearance of isolated surface V O 's at the 2-fold coordinated oxygen (O 2c ) sites.Eventually, the (1 × 2) stripes shown in Fig. 1(a) appear; their density increases upon further reduction until they cover the whole surface [28].Here, we have only isolated (1 × 2) stripes [not a full (1 × 2) superstructure, which would be detectable by diffraction techniques].
This phase separation into (1 × 1) and (1 × 2) regions allows an accurate determination of the critical concentration of V O 's where the transition occurs.In the samples considered here, it is measured as a 16.7 AE 0.2% (≈1=6) monolayer.
A zoom-in of the (1 × 1) phase and a detail of the (1 × 2) reconstruction are shown in Figs.1(f)-1(h) and Figs.1(i)-1(k), respectively.Using an nc-AFM with an O tip termination allows us to image the O sublattice.On the well-known (1 × 1) phase [Fig.1(f)], the image shows the O 2c atoms in the regime of repulsive forces; V O 's appear as missing spots in the rows.This imaging mode is particularly valuable on the (1 × 2) reconstruction [Fig.1(i)], where the structural model is still under debate [4,37,38,[47][48][49][50]. AFM allows a precise measurement of the positions of the surface O atoms in the x-y plane, which helps us exclude most of the proposed structures.Our results confirm the Ti 2 O 3 model [Fig.1(c)] initially proposed by Onishi and Iwasawa [4], recently refined by Wang et al. [37] via DFT-based genetic algorithms, and supported by total-reflection high-energy positron diffraction studies [38].Here, the stripes have a Ti 2 O 3 stoichiometry; i.e., their composition can be viewed as Ti 2 O 4 with 50% of V O 's.The atom labeling for the (1 × 1) and (1 × 2) phases used in the discussion is specified graphically in Figs.1(b) and 1(c).
While STM images [Figs.1(g), 1(h), 1(j), and 1(k)] provide supplementary information about the geometric structure, their importance resides in unraveling the electronic structure of the two phases.In the empty-states STM images of the (1 × 1) phase [the imaging mode mostwidely employed in the literature, Fig. 1(g)], the V O 's appear as bright spots between the surface Ti A S0 atoms.More interesting, however, are filled-states STM images that collect the gap states seen in STS [Fig.1(d)]; the image in Fig. 1(h) directly shows the distribution of trapped polarons originating from the V O 's, as is discussed below.The lattice of the (1 × 2) reconstruction is very flexible and subject to surface buckling; see the DFT ground-state structural model shown in Fig. 1(c).Both experiment and DFT find competing structures with different degrees of buckling (discussed in Ref. [36]).When the buckling is small, like in Figs.1(i Therefore, the number of polarons that form in this structure is in accordance with a nominal V O concentration of 50% (i.e., one missing oxygen per unit cell).This value is also in agreement with the electronic structure results shown in Fig. 1(e): The significant jump in the offstoichiometry from the (1 × 1) to (1 × 2) reconstruction, 16.7% to (locally) 50%, implies a high concentration of polarons and causes the broadening of the in-gap peak in the DOS and filled-states STS [51][52][53] of the reconstructed surface [compare Figs.1(d) and 1(e)].
We now turn our attention to the distribution of polarons in the (1 × 1) structure.It is generally agreed that excess electrons originating from impurities (such as V O 's but also other defects, e.g., interstitial Ti atoms [54]) form polarons preferentially located in the subsurface Ti A S1 [001] rows.These can easily hop to surface sites when thermally activated [17,24,[55][56][57][58], while a polaron trapping deeper in the bulk is significantly less favorable energetically, comparable to delocalized electrons [59].We followed the hopping in FPMD calculations and increased the polaron density by considering slabs with increasing c V O from 5.6% to 50.0%.Indeed, we found that at any c V O , the polarons are mostly located in the topmost two layers, S0 and S1, as reported by the histogram bars in Fig. 2.This result does not depend on the thickness of the simulated slab, corroborating the conclusion that in the (1 × 1) phase excess electrons form near-surface polarons [17,24] and avoid the bulk states suggested in recent literature [60].At the smallest concentration, c V O ¼ 5.6%, the polarons are located in the S1 layer 95% of the time and undergo interand intra-Ti  ), it is energetically more favorable to also host polarons in S0 sites rather than increasing their density in the S1 layer [Fig.2(d)].Interestingly, this c V O value in the calculations agrees with the V O concentration of 16.7%, where the (1 × 1) and (1 × 2) phases are found to coexist in our experiment.Thus, the limited capacity of the subsurface (S1) layer to accommodate the polarons appears to be the driving force behind forming the (1 × 2) reconstruction.
In order to explain this characteristic feature of the polaronic distribution as a function of the concentration of oxygen vacancies, we recall here that polarons are negatively charged quasiparticles.Therefore, their Coulomb-like repulsion should increase with decreasing polaron-polaron distance.Indeed, during the FPMD runs, the probability of finding two polarons at adjacent sites is very low at any c V O .This is seen in Figs.3(c)-3(e), where we report the polaron-polaron autocorrelation functions along the Ti A S1 [001] rows.The highest probability is found for values that maximize the distance between polarons in the row, corresponding to three or four lattice constants apart, depending on the V O density [Figs.3(c)-3(e)].This result is in excellent agreement with experiment, in particular, at c V O ¼ 16.7% where a clear (3 × 1) periodicity is observed in the autocorrelation analysis of filled-states STM images [Fig.3(g)].We notice that, at larger V O concentration (c V O ≥ 22.2%), the peak of the calculated correlation function is not shifted to smaller distances, although the polaron density in the cell increases: The (3 × 1) pattern in S1 is preserved [Fig.3(e)], and polarons populate the S0 layer more often [histogram bar of Fig. 2(d)].This can be understood by inspecting the spatial extension of the polaron charge density.While most of the polaron charge has d z 2 symmetry [see inner yellow isosurface in Fig. 3(b)], a non-negligible amount of charge spreads across the adjacent Ti A S0 and Ti A S1 sites [see external isosurface in Fig. 3(b)].This extended cloud gives rise to the STM contrast observed in experiment [Fig.3(f)]-a characteristic dimerlike shape [Fig.3(a)]-and clarifies the correlation between the distribution of polarons in S1 and the STM spots.The 3 × 1 periodicity of polarons along the Ti A S1 row [sketched in Fig. 3(h)] avoids the overlap of two adjacent extended repulsive polaron clouds [Fig.3(b)], thereby reducing the polaron-polaron repulsion, as recently confirmed in the literature [60].In the FPMD runs at higher V O concentrations, polarons are pushed into S0 sites in order to preserve the optimal concentration in the S1 layer.We also note that the polaron-polaron repulsion is perturbed by the attracting Coulomb interaction between polarons and the immobile V O centers [16,61], which modifies the optimum patterns, including the 3 × 1 (see Ref. [36] for more details).
The measured maximum V O concentration of 16.7% on the (1 × 1) patches of the partially reconstructed surface suggests that a further increase in V O results in a thermodynamically unstable situation, and this leads the system to the (1 × 2) reconstruction.This behavior is borne out of our ab initio calculations, but only if the polaron states are allowed to form [62]: We computed the surface free energy ΔG as a function of the oxygen chemical potential Δμ O for each system [i.e., the (1 × 1) surface with c V O ¼ 5.6%, 11.1%, 16.7%, 22.2%, 33.3%, 38.9%, 50.0% and the (1 × 2) reconstruction].This computation was done both for the unrealistic case of excess electrons delocalized in the conduction band and for the polarons localized in their most favorable configuration.When polaron formation is included in the calculations, the phase diagram in Fig. 4(b) shows that the Ti 2 O 3 reconstruction becomes the most stable one above the critical V O concentration of 22.2%.This is slightly larger than the critical concentration measured experimentally, which could be attributable to effects not included in our calculations, such as the presence of additional dopants (e.g., interstitial Ti [54]), the V O distribution, and the specific type of approximation (DFT þ U) chosen to treat the localization effect [63].Note, however, that the delocalized solution shown in Fig. 4(a) predicts that the reduced (1 × 1) phase remains the most stable structure even at high V O concentrations: The Ti 2 O 3 reconstruction [inset in Fig. 4(a)] is just as stable as the (1 × 1) at c V O ¼ 50.0%.

IV. DISCUSSION
To understand this polaron-mediated change of the surface stability, it is instructive to analyze the energy gained upon polaron formation, E POL [46].It is defined as the total-energy difference between the polaronic and fully delocalized free-carrier solution, and it results from the competition between the strain energy required to distort the lattice, E ST , and the electronic energy E EL gained by localizing the electron at a Ti site in such a distorted lattice.Figure 5(a) shows E POL for the most stable configuration at each considered c V O .For low concentrations, up to c V O ≃ 20% in the (1 × 1) structure, E POL decreases with increasing concentration; it becomes more favorable to form polarons in the system.For c V O > 22.2%, however, E POL remains constant and ultimately increases.Figure 5(a) also includes the polaron energy for the reconstructed surface [64].It has the highest negative value of all configurations considered, and the resulting convex hull clearly shows that the reduced (1 × 1) phases are unstable for c V O > 22.2%.The surface undergoes a structural reconstruction above this critical value.
The trend of E POL as a function of c V O becomes clear if one considers the energy cost that is required to locally distort the lattice in order to accommodate an electron and the energy gained by localizing an electron in such a lattice [46].[001] row.Once the lattice in the S1 layer is distorted (as well as the S0-S1 interlayer distance), it is increasingly easier to form more polarons.The S0 layer is kept essentially unperturbed.At c V O ¼ 22.2%, however, the polarons also start to populate the S0 layer, and the breaking of the symmetry of the top layer results in an additional energy cost.As pointed out above, increased jumping to S0 happens in order to preserve the optimum 3 × 1 pattern in S1 (reached at c V O ¼ 16.7%), as a larger polaron density in S1 would result in the overlap of the extended polaronic cloud, as well as in a strong Coulomb repulsion.While hosting polarons in both layers S1 and S0 reduces this overlap, it requires a distortion of the lattice around both the Ti A S1 and Ti A S0 sites.Apparently, a structural switch to a local (1 × 2) structure happens instead.The Ti 2 O 3 reconstruction exhibits a high degree of structural flexibility, expressed in terms of a low energy cost E ST [Fig.5(b)] needed to accommodate a large concentration of polarons (i.e., as large as in the c V O ¼ 50% case).
Finally, in Fig. 5(c), we show the evolution of E EL as a function of c V O .First, we note that E EL is about twice as larger as E ST , and it leads to a large polaron energy at any V O concentration.To understand the trend of E EL , it is useful to correlate E EL with the change in the polaron DOS with c V O [Figs.1(d) and 1(e), and Fig. S9 in Ref. [36]).At low c V O , the polarons form a well-localized peak below the Fermi energy [Fig.1(d)].At c V O ¼ 22.2%, E EL decreases by 100 meV because the additional polarons start to populate levels at lower energies (Fig. S9).For larger c V O , the peak exhibits a progressive broadening towards higher energies, which leads to an increase of E EL .
The combined action of both the low strain cost E ST and the favorable electronic localization energy E EL determines the thermodynamic stability of the (1 × 2) phase, which overtakes the electrostatic instability from the trapped charges near the (1 × 1) surface.Moreover, the Ti α;β S0 sites in the (1 × 2) reconstructed surface are easily reducible and can conveniently host polarons [see the advantageous formation energy E POL in Fig. 5(a)].The (1 × 2) phase is more stable than the highly reduced (c POL −E 1×2 POL ), and this energy gain arises predominantly from the difference in E EL between the two phases, 64 meV [ΔE ST ð50%Þ ¼ 20 meV].

A. Summary and conclusion
In conclusion, the present study shows that polarons (charges trapped at Ti lattice sites) play an important role for the stability of titania surfaces.At slightly reducing conditions, Ti 3þ polarons are trapped in the subsurface sites of the (1 × 1) TiO 2 surface, with the optimal polaronpolaron distance being three lattice sites.As the polaron density increases, Ti 3þ polarons are pushed to the less favorable Ti surface sites as well.This process, however, costs energy, so ultimately a transition to a (1 × 2) reconstruction with the formal stoichiometry of Ti 2 O 3 occurs.If we neglect charge trapping and relaxation around the polaron, the reconstructed phase does not appear in the theoretical phase diagram.
Our results for the rutile TiO 2 (1 × 1) to (1 × 2) transition could present a new paradigm for surface reconstructions that involve trapped charges and the interaction among them.Polaron formation is ubiquitous in transition-metal oxides and is typically propelled oxygen vacancies and doping [58].Possible examples of other materials in which polaron-mediated reconstructions could be operational are oxide perovskites, in particular, SrTiO 3 [65][66][67], as well as popular oxides like ceria and hafnia [68].Thus, this polaron-mediated mechanism is likely to be a pervasive phenomenon that could explain structural, electronic, and magnetic reconstructions at surfaces and interfaces [69] and could be employed to tune surface properties and to control the surface geometry.For instance, it is known that Nb-doped, V-doped, or Cr-doped TiO 2 samples form the (1 × 2) reconstruction much easier, even when the surrounding (1 × 1) phase has much lower V O concentrations [70][71][72].Finally, the control of charge trapping could provide a way to optimize the functionality of TiO 2 -based memristors [73] and to facilitate charge transfer in catalytic processes [74].
)-1(k), the STM signal comes from the O 3c [empty states, Fig.1(j)] and Ti α;β SO [in-gap states, Fig.1(k)] atoms, the latter is attributed to Ti 3þ polaroniclike states originating from the deviation of the Ti 2 O 3 reconstruction from the stoichiometric formula.Our DFT calculations show that the reconstructed Ti 2 O 3 (1 × 2) surface hosts two polarons per (1 × 2) unit cell.

FIG. 1 .
FIG. 1.The (1 × 1) and (1 × 2) surface phases of TiO 2 ð110Þ.(a) Large-area STM image of a reduced rutile (110) single crystal.Unreconstructed (1 × 1) regions with single oxygen vacancies V O appear together with reconstructed (1 × 2) Ti 2 O 3 stripes.(b,c) Structural models of the (1 × 1) and (1 × 2) phases, respectively.(d,e) STS spectra (line) taken above (1 × 1) and (1 × 2), together with calculated polaronic gap states (filled area), shifted upwards by ≈0.25 eV to align the center of mass of the Ti 3þ (polaronic) peaks with the experiment (in this way, the DFT peak positions are also in better agreement with the calculated vertical excitation energy E EL [46]).(f)-(k) High-resolution images of the (1 × 1) and (1 × 2) regions marked in panel (a), overlaid with the respective structural models.AFM images (i,f) show the highest O atoms (O 2c and O rec ) on the surface, and STM images show the empty (g,j) and filled (h,k) states.
A S1 row hopping [histogram bar of Fig. 2(a)].Hopping into surface sites is very rare, with polarons being trapped in S0 for 5% of the time only.With increasing c V O in the surface (S0) layer [histogram bars of Figs.2(b)-2(d)], i.e., with an increasing number of polarons in the system, hopping into S0 sites becomes more likely.

Figure 2
Figure 2 also presents the charge density isosurfaces of polarons in their energetically most favorable configuration for each c V O .The stable configurations at low V O concentration (c V O ≤ 16.7%) are characterized by polarons trapped only in S1 sites [Figs.2(a)-2(c)].At higher concentrations (c V O ≥ 22.2%), it is energetically more favorable to also host polarons in S0 sites rather than increasing their density in the S1 layer [Fig.2(d)].Interestingly, this c V O value in the calculations agrees with the V O concentration of 16.7%, where the (1 × 1) and (1 × 2) phases are found to coexist in our experiment.Thus, the limited capacity of the subsurface (S1) layer to accommodate the polarons appears to be the driving force behind forming the (1 × 2) reconstruction.In order to explain this characteristic feature of the polaronic distribution as a function of the concentration of oxygen vacancies, we recall here that polarons are negatively charged quasiparticles.Therefore, their Coulomb-like repulsion should increase with decreasing polaron-polaron distance.Indeed, during the FPMD runs, the probability of finding two polarons at adjacent sites is very low at any c V O .This is seen in Figs.3(c)-3(e), where we report the polaron-polaron autocorrelation functions along the Ti A S1 [001] rows.The highest probability is found for values that maximize the distance between polarons in the row, corresponding to three or four lattice constants apart, depending on the V O density [Figs.3(c)-3(e)].This result is in excellent agreement with experiment, in particular, at c V O ¼ 16.7% where a clear (3 × 1) periodicity is observed in the autocorrelation analysis of filled-states STM images [Fig.3(g)].We notice that, at larger V O concentration (c V O ≥ 22.2%), the peak of the calculated correlation function is not shifted to smaller distances, although the polaron density in the cell increases: The (3 × 1) pattern in S1 is preserved [Fig.3(e)], and polarons populate the S0 layer more often [histogram bar of Fig.2(d)].This can be understood by inspecting the spatial extension of the polaron charge density.While most of the polaron charge has d z 2 symmetry [see inner yellow isosurface in Fig.3(b)], a non-negligible amount of charge spreads across the adjacent Ti A S0 and Ti A S1 sites [see external isosurface in Fig.3(b)].This extended cloud gives rise to the STM contrast observed in experiment [Fig.3(f)]-acharacteristic dimerlike shape [Fig.3(a)]-andclarifies the correlation between the distribution of polarons in S1 and the STM spots.The 3 × 1 periodicity of polarons along the Ti A S1 row [sketched in Fig.3(h)] avoids the overlap of two adjacent extended repulsive polaron clouds [Fig.3(b)],thereby reducing the polaron-polaron repulsion, as recently confirmed in the literature[60].In the FPMD runs at higher V O concentrations, polarons are pushed into S0 sites in order to preserve the optimal concentration in the S1 layer.We also note that the polaron-polaron repulsion is perturbed by the attracting Coulomb interaction between polarons and the immobile V O centers[16,61], which

FIG. 2 .
FIG. 2. Polaron distribution in TiO 2 ð110Þ.(a)-(d), Most stable polaron configurations obtained by FPMD and subsequent relaxations at the various V O concentrations specified in each panel.The histograms illustrate how often the polarons are found in S0 (orange) and S1 (green) Ti sites during the FPMD runs.

Figure 5 (
b) reports this strain energy E ST for the most stable configuration at each consideredc V O .E ST decreases with increasing c V O , but it shows a distinct jump at the critical c V O ¼ 22.2% [ΔE ST ¼ E ST ð22.2%Þ − E ST ð16.7%Þ ¼ 65 meV].The isosurfaces plotted in Fig.2are helpful to understand this trend.At low concentrations, up to c V O ¼ 16.7%, the polarons reside preferentially along the Ti A S1

FIG. 3 .
FIG. 3. Polaron distribution along the Ti A S1 row.(a) Simulated filled-states STM (at 2 Å from the surface) of an isolated polaron around a Ti A S1 site.(b) The corresponding charge density isosurface shows the high-density d z 2 -like orbital localized at Ti A S1 (yellow) and the large, low-density cloud extending across the nearby Ti A S1 and Ti A S0 sites.(c)-(e) Calculated polaron-polaron autocorrelation functions along Ti A S1 and Ti A S0 rows extracted from FPMD runs at c V O ¼ 5.6%, 16.7%, and 22.2%.(f) Experimental filled-states STM image at c V O ¼ 16.7%, and (g) corresponding polaron-polaron autocorrelation function along the [001] rows.(h) Sketch of the 3 × 1 periodicity in the Ti A S1 row.

FIG. 4 .
FIG. 4. TiO 2 ð110Þ surface phase diagram.Surface free energy (ΔG) for the most stable configurations obtained from the MD for the reduced (1 × 1) structures at different V O concentrations (all lines except pink) and for the (1 × 2) Ti 2 O 3 reconstruction (pink line) as a function of the chemical potential Δμ O with respect to the stoichiometric (1 × 1) phase.(a) Delocalized solution: Excess electrons are spread across all lattice sites and delocalized at the bottom of the conduction band.At c V O ¼ 50.0%, the (1 × 1) structure is degenerate with the (1 × 2) Ti 2 O 3 reconstruction (Ti 2 O 4 with 50% of V O ).(b) Localized solution: The excess electrons are trapped in distinct Ti sites, mostly located in S0.The insets show the distribution of excess electrons in the (1 × 2) structure.The localized polaron picture is capable of predicting the (1 × 1) to (1 × 2) transition, manifested by the increased stability of the reconstructed Ti 2 O 3 phase at low oxygen chemical potential compared to the highly reduced (1 × 1) phases.