Machine Learning Based Prediction of Polaron-Vacancy Patterns on the TiO$_2$(110) Surface

The multifaceted physics of oxides is shaped by their composition and the presence of defects, which are often accompanied by the formation of polarons. The simultaneous presence of polarons and defects, and their complex interactions, pose challenges for first-principles simulations and experimental techniques. In this study, we leverage machine learning and a first-principles database to analyze the distribution of surface oxygen vacancies (V$_{\rm O}$) and induced small polarons on rutile TiO$_2$(110), effectively disentangling the interactions between polarons and defects. By combining neural-network supervised learning and simulated annealing, we elucidate the inhomogeneous V$_{\rm O}$ distribution observed in scanning probe microscopy (SPM). Our innovative approach allows us to understand and predict defective surface patterns at previously inaccessible length scales, identifying the specific role of individual types of defects. Specifically, surface-polaron-stabilizing V$_{\rm O}$-configurations are identified, which could have consequences for surface reactivity.


Introduction
6][7] In the specific case of so-called small polarons, the polaronic charge is localized almost entirely on one atomic site, surrounded by sizable distortion of the local lattice structure. 83][24] The localized charge carriers act as active centers, which enhance (photo)catalytic activity by providing sites that can readily adsorb and interact with reactant molecules. 25,26Although polaron formation may in principle occur on any site of the lattice, the defects can act as attractive or repulsive centers, favoring specific polaronic configurations over others. 27In turn, the dynamics and distribution of the atomic defects are known to be altered by the polarons. 28Therefore, control over the spatial distribution of polaronic active centers becomes pivotal in optimizing (photo)catalytic performance.
While theoretical studies based on density functional theory (DFT) have elucidated excess charge localization in relation to the inducing defect in many materials, [29][30][31] the specific role of subsurface and surface polarons, particularly in the presence of defects, on the archetypal redox active oxide surface TiO 2 (110) is still debated.Here, a problem arises from the complexity of the configuration space of point impurities, where DFT calculations strive to account for the computational cost of the problem.As a consequence, either no exploration attempt is performed (i.e., most studies rely on the configuration randomly obtained in the DFT calculation), 32,33 or effective but costly approaches are adopted such as molecular dynamics, 34,35 Monte-Carlo-driven DFT simulations 36 or systematic explorations limited to a handful of localization sites. 37Thus, finding a method that effectively navigates the diverse defect-polaron configuration landscape has become a research imperative.
In this study, we focus on rutile TiO 2 (110) and show how the spatial distribution of V O measured by SPM can be successfully predicted and interpreted by first-principles calculations if the coupling between V O and polarons is taken into account.To address this problem, we developed a strategy based on defect distribution descriptors and neural networks to predict the stability of specific polaron-vacancy patterns.Through an iterative optimization active learning cycle, we systematically extended the DFT reference dataset and converged the machine learning (ML) model, to efficiently explore the defect-polaron configuration space.The model can capture the complexity of the V O -polaron interactions with DFT accuracy and proposes new configurations showing remarkable energy stability.
By feeding Markov-chain Monte-Carlo (MC) algorithms with the ML configuration energies, we simulate the annealing process leading to the formation of vacancies and polarons in the experimental samples.As a final result, we obtain large-area (>10×10 nm 2 ) surface morphologies resembling the SPM measurements.This analysis revealed new physical properties of the polarons on TiO 2 (110), where the formation of inhomogeneously distributed V O is linked to an increased formation of surface polarons and, therefore, to the density of active sites.

Results and Discussion
Defect Distribution via DFT, Experiment, and Machine Learning Fig. 1 shows the surface structure of reduced rutile TiO 2 (110) as imaged by constant current STM measurements (see panel b and Methods Section), together with the models predicted from DFT without taking polarons into account or by explicitly modeling their impact via machine learning (see panels a and c respectively).The unreconstructed 1×1 rutile surface consists of alternating rows of under-coordinated (two-fold) oxygen atoms (the bridging oxygen atoms, O br ) and five-fold coordinated titanium atoms (Ti 5c ) running along the [001] direction. 38,39Oxygen vacancies form easily on the O br sites upon sputtering and annealing, up to a critical concentration of c V O ≃ 17%. 341][42][43][44][45] Every V O releases two excess electrons that form polaronic states, localizing preferably on subsurface Ti sites. 27,35,46Thus, the V O can be considered as a positively charged (2+) center.By simple electrostatic considerations (and by, simultaneously, neglecting the role of polarons), one would expect a purely repulsive interaction among the vacancies.In this picture, the configuration maximizing the V O -V O distances represents the most favorable vacancy distribution.For the critical concentration of c V O = 17%, this corresponds to a homogeneous configuration with a V O -V O distance of 6 lattice sites along the [001] row, and 3 lattice sites considering two oxygen vacancies on adjacent rows (see Fig. 1a).
DFT calculations confirm the homogeneous V O distribution in Fig. 1a as the ground state configuration, as far as the formation of the polarons is suppressed (i.e., the excess electrons are forced into spatially delocalized states at the bottom of the conduction band, rather than localized polaronic states).While this unphysical metallic solution (rutile TiO 2 is an n-type semiconductor) is less stable than the polaronic solution, it simplifies the search for the optimal defect distribution via a two-step process.Initially identifying the optimal defect pattern through DFT calculations, where polaron formation is suppressed, and subsequently introducing polarons into random positions or finding the most favorable polaron configuration within the given defect distribution. 27While this approach reduces the combinatorial divergence of defect-polaron configurations, it relies on the assumption that the distribution of atomic defects is not affected by the polarons, which is not valid for most materials. 28e experimental measurements do not support such homogeneous V O pattern.Fig. 1b shows a typical image as obtained from low-temperature STM measurements on a TiO 2 (110) surface after sputtering and annealing treatment to form a high content of oxygen vacancies (c V O ≃ 14%, close to the critical value of 17%).At this temperature, the oxygen vacancies (imaged as bright spots along the dark [001] O br rows) are immobile and appear in irregular patterns, quite far from any homogeneous distribution.The discrepancy with the simple models discussed above is a strong indication of the role that polarons can have in determining the optimal V O surface structure.Simply adding polarons on a rigid V O pattern (effectively decoupling V O and polaron) as usually done in standard DFT simulation, would not improve the situation.Fig. 1c reports the surface structure as predicted by our machine learning model, which allows simultaneously varying both V O and polaron positions to find the configuration that minimizes the total energy of the system.The resulting V O distribution is in good qualitative agreement with the inhomogeneous distribution found throughout experiments.Our methodology, described in detail in the following, is capable of capturing the effects of the polarons on the oxygen vacancy distribution, going beyond the simple picture relying on purely V O -V O interactions.Moreover, it allows us to consider large surface areas of about 250 nm 2 (>15×15 nm 2 ), corresponding to 54×24 supercells, extending considerably the limits of standard DFT simulations.

Machine Learning Polaron and Defect Distributions
The methodology proposed here is structured in three parts: First, we train a feed-forward neural network 47 to predict the DFT energy of the system depending on the configurations of the impurities.Due to the computational limitation of the DFT calculations, we adopt relatively small unit cells in this step.Specifically, we used two supercells with different lateral extensions (6×4 and 12×2, see Methods Section) to include long-range interactions along different crystallographic directions.Then, we use the trained model to search for low-energy configurations that were not included in the original set of data, adopting an active learning scheme. 48Finally, we use the actively trained model to obtain large-area predictions.In the following, we describe in detail the architecture of the machine-learning model and compare the ML predictions with experimental data on reduced TiO 2 (110).
The training of the machine learning model requires a reference database built up by several, distinct polaron and atomic-defect configurations.By following the process described in detail in the Methods Section, we calculated the free energy for different configurations at the DFT+U level using VASP, [49][50][51] with a U = 3.9 eV on the d orbitals of Ti atoms. 27,52larons were localized at chosen surface Ti S0 (Pol S0 ) and subsurface Ti S1 (Pol S1 ) sites via occupation matrix control. 53We modeled 2367 symmetrically-inequivalent polaron-V O configurations in a 6×4 supercell (i.e., six and four times the  surface is depicted in Fig. 2a.The descriptor representing the configuration is constructed by, first, discretizing the space into a rough grid (see Fig. 2b), and encoding the spatial distribution of polarons and atomic defects.The discretized space simplifies the training of the ML model, as compared to using Euclidean distances. 27To improve the description of spatially close polarons/defects, we employ a one-hot encoding (i.e., value of 1 for grid cells containing a defect, 0 otherwise), smeared via multiple applications of a discrete Laplacian kernel (Fig. 2c).Then, to predict the energy of the whole system given a specific configuration, we split the total energy into contributions arising from a single defect/polaron impurity (Fig. 2d): Here, E tot is the total energy of a given configuration, and E i , E j and E k are the virtual contributions of a single Pol S1 , Pol S0 , and V O respectively.We use a feed-forward neural network to estimate the virtual contribution of a single defect/polaron (Fig. 2e).Finally, we sum over the virtual contributions to obtain the total configuration energy. 27The total energy can be computed by DFT calculations, 37 while the virtual contributions are not directly accessible in the DFT data.Thus, we can train our model using the discretized Aiming for a comparison with the experimental measurements, we focus here on the lowenergy configurations, which are more likely to get stabilized in real samples.To identify such stable configurations, we performed simulations that model the annealing process.In the preparation of the experimental samples, both polarons and oxygen vacancies diffuse on the sample during annealing.At lower temperatures, V O on rutile are immobile, while polarons always show a certain degree of mobility, hopping/tunneling a few lattice sites around the equilibrium position. 34,39The simulated annealing can be implemented as a global optimization scheme. 54Candidate configurations are obtained by perturbing the current configuration, randomly displacing one defect/polaron impurity to any nearest neighbor site.
The new configuration is either accepted or declined by virtue of the Metropolis-Hastings algorithm 55 with the acceptance criterion based on the configuration energy -similar in spirit to large-scale defect distribution studies based on reverse MC 56 (although in our approach the defect distribution is not fitted to minimize the deviation from experiment, but it relies entirely on DFT/ML data).Calculating the energy of the candidate configurations within the DFT framework would make this approach unfeasible, due to the computational cost of DFT calculations and the high number of energy evaluations required for a single optimization.
Conversely, the ML model allows us to inspect the stability of an extremely high number of defect-polaron configurations (minimization of the energy requires on the order of 10 3 − 10 6 energy evaluations depending on the size and initialization of the employed configuration) and enables the exploration of candidate structures.
We iterated annealing simulations following an active learning procedure.The initial DFT data set built by random configurations was progressively augmented by including the results from the annealing optimization (see SM Fig. S3).Consequently, we obtained a final ML model refined to account for a broader range of configurations.The refined model is finally used to obtain large-area predictions (7.1×10.5 nm 2 to collect statistics and 16×16 nm 2 for visualizations) on the defect-polaron distributions, using again the simulated annealing approach.The qualitative agreement with the experimental data is shown in Fig. 1c for the 54×24 rutile TiO 2 (110) supercell.In the following, we quantitatively analyze our results.

Formation of V O -Polaron Patterns and Their Mutual Interaction
The analysis of the low-energy configurations (see SM Fig. S4 for the energy distribution of all possible V O configurations in the 6×4 cell) is summarized in Fig. 3. Fig. 3a shows the improvement of energies of the TiO Full inclusion of polaron-V O interaction via our proposed ML protocol ("ML Polarons").
By suppressing polaron formation, the ground state configuration is given by the vacancies being homogeneously distributed on the surface (C Hom NoP configuration, see "No Polarons" column in Fig. 3a).The "Random Polarons" column of Fig. 3a shows instead the energy of the system obtained by including polarons in random positions and enriched by adding specific, low-energy polaronic configurations that were suggested in previous studies. 27,34,37re, the ground state configuration is given by a homogeneous distribution of Pol S1 in the homogeneous V O pattern (labeled as C Hom RandP in Fig. 3a).Treating polaron-V O coupling at the ML level (third approach) results in novel V O distributions with lower energy, indicating a new ground state for the system, where the homogeneous configuration is no longer the most stable one, as shown in the "ML Polarons" column of Fig. 3a.First, we note that the ML model identified a different order of Pol S1 showing a better stability in the homogeneous V O background (labeled as C Hom ML , see also SM Fig. S5).Moreover, new polaron configurations explored by the extensive ML search improve the stability of many other V O patterns (see the energy levels in black in the "ML Polarons" column of Fig. 3a, lower than in the "Random Polarons" column).Importantly, two of these previously-unexplored polaron configurations (labeled as C 0 ML and C 2 ML ) resulted in energy values even lower than the homogeneous distribution, revealing a new ground state for the system.Moreover, novel V O -patterns were proposed by the ML search as low-energy configurations.One in particular (red line in Fig. 3a) is ranked as the second most stable configuration (C 1 ML ).The polarons play a key role in stabilizing this V O -pattern and as further proof, we calculated the energy of this new V O -pattern, artificially suppressing the polaron formation, and obtained a much worse stability (red line in the "No Polarons" column).
Interestingly, in all the new low-energy configurations obtained from the ML-driven search (except for C Hom ML ), we note the presence of at least one polaron on a surface Ti S0 site (configurations containing Pol S0 are orange highlighted in Fig. 3a).Fig. 3b and c compare the spatial distribution of the surface Pol S0 and subsurface Pol S1 .The formation of the surface polaron is particularly stable when occurring in the central Ti S0 site between two oxygen vacancies aligned on the [ 1 12] direction (see top view in Fig. 3b).This [1 12]-aligned V O -polaron complex represents indeed the ground state configuration obtained by our ML search (e.g., it is present in C 0,1,2 ML ).Another remarkably stable complex is given by two vacancies aligned along the [1 10] direction and one Pol S0 in their vicinity (see SM Fig. S5).For instance, this complex appears in the configurations highlighted in purple in Fig. 3a (C 10,11 ML ).The [1 12] and [ 1 10] alignments found in the ML search agree well with the experimental SPM measurements (compare Fig. 1b and c) showing a high coverage of such high density V O -regions.
In contrast, DFT predictions, which neglect polaron-V O interaction or randomly distribute polarons, favor homogeneous configurations.

Comparison of a Large Scale Model and the Experimental Surface
Fig. 4 shows our results as obtained by ML-driven annealing simulations on large-area 24×16 cells (corresponding to 7×10 nm 2 ), which enables a direct comparison with the experiment.
Visual inspection (Fig. 1b vs. c) already indicates that our ML treatment provides V O distribution that closely resembles the experimental one.We quantify this agreement by calculating autocorrelation functions (ACF) 57 for simulated annealings under different computational conditions and compare it to the experimental ACF of the V O distributions extracted from Fig. 1b (for details see Supplementary Fig. S7).The simulated annealing procedure starts from random V O -polaron configurations, where we obtain several large-area models (such as the one in Fig. 1c), all showing very similar characteristics.To complete our comparison, we also use the ML model to anneal a system where polaron formation is suppressed.This scheme, similar to the non-polaronic DFT approach of Fig. 1a, assumes a homogeneous V O pattern but takes into account annealing-induced disorder effects.
The ACFs are shown in Fig. 4, where projections of V O defect populations along the same and adjacent [001] rows are shown in the histograms in panels a and b, respectively.
For oxygen vacancies lying on the same row, both the ML model and the experiments show that short V O -V O distances of 1 and 2 lattice sites are unlikely.The highest probability lies at a distance of 4 or 5 lattice sites for both the experimental and ML annealing including polaron-V O interactions (see Fig. 4a and b red and blue data, respectively).By considering only the V O -V O repulsion as driving force (i.e., excluding polaron formation in the ML annealing procedure; see ML No Pol in Fig. 4a and b) and applying an identical annealing protocol as in the polaron-V O interaction case, we find the probability maximum lying at a 6-site distance for in-row and 3-site distance in the adjacent row.This is further evidence for polarons' role in stabilizing the V O arrangement.
As a result, the rutile TiO 2 surface shows some areas with a locally low density of oxygen vacancies (down to 0%), alternated with highly dense areas (up to 20%, which is compatible with the 4-site-distance distribution).Our data suggest that the great stability of the [1 12]and [1 10]-aligned V O -polaron complexes contribute to this alternation of locally less and more reduced areas at a given c V O .To further corroborate this result, we performed additional DFT calculations modeling this strong inhomogeneity (see SM Fig. S5 with configurations C 26−29 ML ).We also note that this analysis reconciles the DFT predictions on the critical concentration at which the (1×2) surface reconstruction occurs for the surface phase transition, which was calculated as ∼ 20%, in apparent disagreement with the experiments reporting an average concentration of 17%. 34

Conclusions
In summary, we directly elucidated the impact of polarons on the structure of oxide surfaces, using an example of the prototypical rutile TiO 2 (110) surface.Specifically, we designed a computational machinery to predict the distribution of polarons and oxygen vacancies on rutile TiO 2 (110), by performing machine-learning-guided DFT calculations.MC-driven annealing simulations based on the ML data enabled the exploration of defect distributions on scales much larger than standard DFT allows.An analysis of the experimental SPM images yielded a direct validation of the theoretical predictions.While conventional approximations used in traditional DFT calculations result in homogeneous solutions, we were able to retrieve the inhomogeneity of the V O distribution as detected by the experiments.Our analysis clarifies the peculiar inhomogeneous distribution of V O on rutile TiO 2 (110).Most importantly, the system shows a tendency towards the formation of high-density V O patterns, alternated with low-density V O regions.While larger defect-free areas are typically attributed to subsurface Ar impurities, 58 the here observed fluctuation of the local c V O can partially be attributed to the interaction of polarons and V O s.
These results suggest that surface reactivity could be optimized by tuning the annealing procedure to facilitate the formation of energetically more favorable, high-density V O patterns, which promote surface localized charges and their interaction with adsorbates. 26To elucidate the role of the surface polaron, further experiments are necessary.Resonant photoelectron diffraction does not rule out the formation of surface localized charge carriers, even at low c V O . 59SPM measurements in the presence of CO adsorbates confirm the formation of the [1 12]-aligned V O -polaron complex, 26 while STM measurements probing the filled states on the clean surface do show some disparity in comparison to simulated STM. 37The reasons for this discrepancy are manifold, ranging from temperature-induced effects, 60 to the electric field of the tip.
We expect our methodology to be applicable to any other polaronic system, even including multiple defects as sources of polarons, such as the perovskite SrTiO 3 (001) surface 61 exhibiting Sr adatom/vacancy and often doped by Nb atoms. 29Moreover, this methodology could be used to study the spatial distribution of defects (e.g., subsurface, bulk) that are not directly accessible by the experiments, such as interstitial titanium in rutile.Additionally, the stochastic optimization model could be further improved by considering realistic anisotropic diffusion probabilities along certain directions.This could be achieved by explicitly computing hopping and diffusion barriers, and incorporating these barriers into the annealing simulations.

DFT Modeling
We performed DFT+U calculations using VASP [49][50][51] on the rutile TiO 2 (110) surface.We used standard projector augmented wave pseudopotentials for Ti (treating d-and s-orbitals as valence) and soft Oxygen pseudopotentials.We adopted a Hubbard U = 3.9 eV on the d orbitals of Ti atoms. 27,52The sampling of the reciprocal space included the Γ-point and the plane-wave energy cutoff was set to 400 eV.
The surfaces were modeled using 5-layer-thick slabs (where the two bottom stochiometric layers were fixed at their bulk position) with lateral supercell sizes of 6×4 and 12×2.
To partially account for the role of thermal effects in the stabilization of the V O patterns during the annealing treatment in the experiments, we used an expanded [001] lattice vector.Specifically, the low T lattice constant of 2.953 Å38 was expanded to 2.968 Å (high T corresponding to 500-600K) in accordance with thermal expansion coefficient measurements. 62,63This strain of +0.5% is well below the crossover point of +3%, where surface polaron formation is favored over subsurface polaron formation. 60thin the supercells, we removed 4 surface bridging oxygen atoms (in random positions) from every slab, obtaining a c V O of approximately 17%.To assess non-polaronic solutions we performed spin unpolarized DFT, constraining the excess electrons in spatially delocal-ized states at the bottom of the conduction band.To model the polaronic structure, we followed a three-step procedure: Initially, we removed bridging oxygen atoms from a pristine structure to generate a specific oxygen vacancy pattern.This structure was relaxed while all excess charge carriers were kept delocalized by employing a spin un-polarized relaxation.
After retrieving the structural properties of the oxygen vacancy configuration, we introduced polaronic distortions at selected sites via occupation matrix control, 53 using distinct occupation matrices for Pol S1 and Pol S0 sites. 27Finally, we performed an unconstrained relaxation starting from the structures and wave functions determined in the previous step.

ML Model Training and Optimization of Defect Configurations
The machine learning model is implemented in the framework of JAX. 64Here, we describe the model optimization based on the study of configurations in the 6×4 supercell.We optimized the machine learning model using stochastic gradient descent and backpropagation on an augmented dataset, including all symmetrically equivalent representations, of the training defect configurations.We randomly split this dataset into 80% training data and 20% validation data and optimized the model parameters by minimizing the mean squared error of the energy prediction of the training data via backpropagation.Before training, energies, as well as the descriptors, were rescaled to [0, 1], by min-max scaling.Using an early stopping mechanism, the best model was selected based on the lowest validation dataset error within the optimization procedure.The mean squared error during training as well as a scatter plot of DFT and ML energies are displayed in Supplementary Fig. S1 and S2.
To ensure sufficient accuracy when using the model in the case of exploration, we applied an active learning procedure as depicted in Supplementary Fig. S3.Here, we performed an iterative training-testing loop to further improve the reliability, data efficiency, and scope of the proposed model.Since our main interest lies in the determination of low-energy polaron-defect complexes, our model was used for the optimization of defect configurations in various cases.We searched for global optima of configurations by allowing all defects to diffuse during the optimization.Local minima of fixed polaron layer densities were added by restricting polaron movement to intra-layer hopping.Also, local minima of cases where the V O -configuration was fixed and only polarons were relaxed, were explored.Within these three exploration cases, we extracted and confirmed the most stable configurations by performing a comparative DFT calculation of the proposed polaron configurations.
The optimization of configurations is performed via simulated annealing, 54 where the temperature variable in the Metropolis criterion was set to 1000K (similar to the annealing temperatures in the sample preparation).Even though the diffusion processes of the respective defects during the optimization are physically motivated, they do not necessarily represent the physical process of the formation of observed defect patterns.Defect transport mechanisms such as inter-row hopping of oxygen vacancies have not been reported 65 but may improve optimization efficiency or more efficiently overcome energy barriers.Discrepancies between polaron and vacancy hopping rates were also ignored, which potentially affects the final outcome of the optimization.Similar effects were observed for the specific temperature or temperature ramp employed in the simulated annealing.

Experimental Setup
SPM was performed using STM in an ultrahigh vacuum (UHV) chamber with a base pressure below 2×10 −11 mbar; the whole chamber, equipped with an Omicron qPlus low-temperature head, was suspended using 36 bungee cords for efficient vibration damping. 66 Figure 1b displays the z-channel of a feedback-controlled unoccupied-states STM image taken at a sample temperature of 14 K; Imaging parameters: sample bias V S = +0.9V, grounded tip, tunneling current set-point I t =20 pA, oscillation amplitude A=500 pm.
The contrast in Fig. 1b corresponds to a typical unoccupied-states STM imaging contrast over a reduced rutile TiO 2 (110) surface, which is dominated by electronic rather than geometric considerations: 1 eV above the Fermi level the conduction band consists of Ti 3d states and defect V O states, while the O states constitute the valence band. 70Therefore, the highest probability of electron tunneling from the tip to the surface is above the Ti 5c rows and V O s -they appear bright under these STM conditions.On the other hand, the tunneling is less likely above the O br rows and they appear dark even though they geometrically protrude highest from the surface.Note that in Fig. 1b individual Ti 5c atoms can be recognized as spheres forming a row along the [001] direction, while V O s are recognized as isolated, bright spheres.During the active learning scheme, we exhaustively predicted optimal polaron configurations in all possible oxygen vacancy configurations in the 6 × 4-cell.Figure S4 shows the energy distribution of the most favorable configurations as predicted by the ML model.The marked region at the lowest energies shows the most interesting configurations and is comprised of all configurations that feature one oxygen vacancy per [001] bridging oxygen row.Less stable configurations were not explored as thoroughly during active learning and might contain less accurate predictions, or some configurations, which were not reproducible at the DFT level.
In Figure S5, we collect the best polaron configuration for all V O -configurations in the low energy cluster, with their respective DFT energy.Figure S6    cluster.Some notable features here are the presence of Ti S0 -polarons (marked in orange) as compared to the reference dataset, where most stable configurations usually only consisted of Ti S1 -polarons (marked in yellow).The presence of Ti S0 -polarons is usually associated with two oxygen vacancies aligned along [1 12] or [ 1 10], which stabilizes the excess charge at a surface Ti 5c site.Subsurface polarons are most stable in the homogeneous ground state, in particular, when placed symmetrically and diagonally aligned between two oxygen vacancies in adjacent rows.First, we identified the V O positions in the STM image (Figure S7a).Then, we calculated the average number of V O s around each detected V O position.This was done by overlaying copies of the V O positions on top of each other, such that each V O was placed at the origin once (Figure S7b).The resulting grid-like distribution was subdivided into individual grid cells, and the number of V O s in each cell was counted.

Figure 1 :
Figure 1: Oxygen vacancy distribution on rutile TiO 2 (110) obtained by various methods.(a) Schematic representation of the most favorable V O distribution in non-polaronic DFT calculations as obtained from a 6×4 (∼1.8×2.6 nm 2 ) supercell.The schematic depiction is generated by showing the O br bridging atoms as black regions and Ti 5c rows and V O as white.The inset displays the structural model of rutile TiO 2 (110).The distance maximizing V O distribution (6 sites in row, 3 sites in adjacent row) and the 6×4 supercell are indicated.(b) Unoccupied-states, constant-current STM image of a clean, reduced rutile TiO 2 (110) surface (imaging parameters in the Figure) depicting Ti 5c rows and V O s as bright, while O br rows are depicted as dark.More details on the contrast formation are given in the Methods.Locally low and high V O -concentrations (c V O ) areas are marked with solid and dashed red boxes, respectively.The crystalographic directions are consistent in all panels.(c) ML-predicted schematic representation of surface oxygen vacancy distribution, where the interaction of surface and subsurface polarons (Pol S0 and Pol S1 , respectively) and V O s are modeled in a 54×24 (∼ 16×16 nm 2 ) supercell.Orange and yellow markers show the position of surface and subsurface polarons in the ML prediction.
figurations in a 6×4 supercell (i.e., six and four times the [001] and [1 10] lattice vectors, respectively), and 2155 configurations in a 12×2 supercell.To optimize the model, we randomly split the calculated configurations and energies into training and validation data sets, including 80% and 20% configurations, respectively.

Figure 2 :
Figure 2: Machine learning model architecture.(a) A defect structure consisting of oxygen vacancies and polarons in a supercell.(b) The supercell is converted into a discretized grid, where each cell encodes whether it contains a defect/polaron.(c) Smearing of the one-hot encoding.(d) The supercell is partitioned into the local environment of each defect.(e) The local environment descriptors are fed through a feed-forward neural network to predict the energy contribution of each defect.The sum of the individual defect contributions gives the total energy of the system.
defect-polaron positions as a descriptor, and the DFT energy as the target quantity.By training the ML model on DFT data obtained for the 6×4 unit cell (see SM Fig. S1), we achieved a mean absolute error of 1.8 and 2.2 meV/V O for the training and validation sets, respectively.By adding training data from the 12×2 unit cell (see SM Fig. S2), the mean absolute error increased slightly (2.9 and 3.5 meV/V O in training and validation, respectively).However, by using both sets of data in the training, the ML model can account for longer interactions in both the [001] and [1 10] directions.For a detailed description of the training process see Methods Section.

Figure 3 :
Figure 3: Analysis of V O -polaron configurations in a 6×4 TiO 2 (110) cell.(a) Comparison of selected low-energy vacancy-polaron configurations as obtained by different treatments of the polaron-V O interaction.For a comparison of all configurations and their labeling, refer to Supplementary Fig. S5.The change in energy for all low-energy configurations is displayed in Supplementary Fig. S6."No Polarons" refers to DFT calculations suppressing polaron formation."Random Polarons" refers to the reference DFT data set, built by including polarons in random positions or guided by physical intuition."ML Polarons" indicates the DFT energies of configurations identified in the ML search.Total energies ∆E are shown using the homogeneous V O distributions (from "No Polarons" and "Random Polarons") as references (note the large energy gain ∆E pol of -3.23 eV between non-polaronic and polaronic solutions with homogeneous V O patterns).Dashed lines connect identical V O configurations.New V O configurations found in the ML search are displayed in red.The occurrence of Ti S0 polarons is highlighted in orange and purple for [1 12]and [1 10]-aligned oxygen vacancies, respectively.The most important V O -polaron complexes are shown schematically in top view at the bottom of the Figure.Only the most stable polaronic configuration per V O arrangement is shown.(b,c) Top and side views of the polaronic isocharge surfaces of the [1 12]-aligned V O -Pol S0 complex (top), and of the Pol S1 in the homogeneous V O -distribution (bottom).

Figure 4 :
Figure 4: Autocorrelation functions of the V O positions as extracted from Fig. 1b (Exp) and ML-based ACFs with (ML Pol) and without polarons (ML No Pol).(a,b) Comparison of experimental and simulated V O autocorrelation functions along a single (a) and adjacent (b) [001]-aligned O br row.Experimental autocorrelation functions are averaged to remove remaining anisotropies.Simulated ones are averaged over all symmetrically equivalent most stable configurations from 60 differently seeded simulated annealing runs performed in 24×16 supercells.ML Pol and ML No Pol are started from identically seeded V O patterns, with 2N V O polarons and no polarons, respectively.Autocorrelation functions are rescaled to account for c V O differences in experiment (14.2%) and simulation (16.7%).
Stiff qPlus sensors67 (k = 1800 N•m −1 , Q=5000-30000, f 0 ∈[25-45] kHz) with a a sharp W tip68 were used to collect the tunneling current (I t ) and the frequency shift (∆f ) signals; deflection detection was achieved using a cryogenic preamplifier in vacuum.69W tips were treated at a Cu(110) surface decorated with a sharp, conductive Cu pyramid at the apex, and were subsequently applied for imaging the rutile TiO 2 (110) surface.Tip sharpness was indicated by the low frequency shifts (∆f ∈ [0, −1] Hz) recorded during STM imaging of a Cu(110) test sample.Sample preparation was performed in a separate UHV chamber (connected to the measurement chamber via a gate valve for in-situ transfer) with a base pressure below 1×10 −10 mbar.Surfaces were cleaned by cycles of sputtering and UHV annealing that consequently reduced the samples and introduced V O s to the surface.A typical cleaning cycle consisted of sputtering with 1.5 keV Ar + ions for 10 min with an ion current of 1 µA• cm −2 , and subsequently annealing the sputtered surfaces in UHV up to 700 • C. Before each measurement, 3-5 cleaning cycles were performed.The over-reduction of the surface was occasionally remedied by annealing the sample to 750 • C in 5×10 −7 mbar of O 2 shower for 10 min.When the reduction level was too high, the rutile TiO 2 samples were re-oxidized ex − situ at 800 • C in O 2 flow and reintroduced to UHV for cleaning.

Figure S3 :
Figure S3: Comparison of test-set predictions for the model with (up-pointing triangles) and without (down-pointing triangles) iterative active learning.Blue and green data points were not added to the training data, orange and red have been used to extend training data.The inset schematically displays the active learning loop Figure S4: Distribution of ML-predicted lowest energies within all possible V O -arrangements in the 6×4-cell.Only energies of the optimized polaronic structure for each V O -configuration are displayed.The highlighted cluster in the bottom left corresponds to the in-depth optimized configurations displayed in Figure S5.

Figure S5 :
Figure S5: Schematic representation of the ML-determined optimal polaron-V O configuration.The 30 most stable V O -configruations are displayed as determined by simulated annealing of the polaron configuration in a fixed V O -configuration.V O are displayed as dashed red circles, surface polarons in orange, and subsurface polarons in yellow.The DFT energies are given relative to the random polaron homogeneous V O distribution.The number in the label indicates the stability rank.

Figure S6 :
Figure S6: Changes of energies for all V O configurations in the low energy cluster by different treatments of polaron-defect interactions.The ML Polarons column corresponds to the configurations displayed in Figure S5, where the most stable configuration in the bottom of the ML Polarons column is c 0 ML and the highest energy level belongs to c 29 ML .The other configurations are labeled sequentially.Configurations containing S0 polarons are highlighted in orange.The left energy scale ∆E deloc labels the No Polarons column.The other two polaronic columns are labeled by the right energy scale ∆E loc .

Figure S7 :
Figure S7: Experimental autocorrelation function as extracted from Fig. 1b.(a) Detected V O positions marked by red circles in STM measurement.(b) Experimental autocorrelation function of the V O positions along a single and in the adjacent [001]-aligned O br rows.