# Machine Learning Based Prediction of Polaron-Vacancy Patterns on the
TiO2(110) Surface

Viktor C. Birschitzky<EMAIL_ADDRESS>[  Igor Sokolović  [
Michael Prezzi  [  Krisztián Palotás  [  Martin Setvín  [  Ulrike Diebold  [
Michele Reticcioli  [  Cesare Franchini  [

###### Abstract

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 (VO) and induced small polarons on rutile TiO2(110), effectively
disentangling the interactions between polarons and defects. By combining
neural-network supervised learning and simulated annealing, we elucidate the
inhomogeneous VO 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
VO-configurations are identified, which could have consequences for surface
reactivity.

UniVie] Faculty of Physics and Center for Computational Materials Science,
University of Vienna, Vienna, Austria TU] Institute of Applied Physics, TU
Wien, Vienna, Austria UniVie] Faculty of Physics and Center for Computational
Materials Science, University of Vienna, Vienna, Austria UniBud] Institute for
Solid State Physics and Optics, HUN-REN Wigner Research Center for Physics,
Budapest, Hungary UniPra] Department of Surface and Plasma Science, Faculty of
Mathematics and Physics, Charles University, Prague, Czech Republic [TU]
Institute of Applied Physics, TU Wien, Vienna, Austria TU] Institute of
Applied Physics, TU Wien, Vienna, Austria UniVie] Faculty of Physics and
Center for Computational Materials Science, University of Vienna, Vienna,
Austria UniVie] Faculty of Physics and Center for Computational Materials
Science, University of Vienna, Vienna, Austria [UniBo] Dipartimento di Fisica
e Astronomia, Università di Bologna, Bologna, Italia

##  1 Introduction

The rich and tunable physics of oxides depend on their precise chemical
composition, and the presence of impurities, including atomic vacancies,
interstitial atoms, and dopants in the material 1, 2, 3, 4. Defects at the
atomic level frequently lead to the formation of polarons, which are localized
charge carriers arising from the synergy between unbound charges and lattice
phonons 5, 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 8. In conjunction with
their inducing defects, these small polarons play a dominant role in a wide
range of processes relevant to technological applications 9 and fundamental
phenomena such as charge carrier mobility 10, 11, 12, electron-hole
recombination13, 14 and adsorption15, 16, 17.

Most importantly the role of polarons is known to be highly relevant in the
context of (photo)catalysis 18, 19, 20, 21 and, single-atom catalysis 22, 23,
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, 26. Although 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
27. In turn, the dynamics and distribution of the atomic defects are known to
be altered by the polarons 28. Therefore, 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 TiO2(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 37. Thus, finding a method that effectively
navigates the diverse defect-polaron configuration landscape has become a
research imperative.

In this study, we focus on rutile TiO2(110) and show how the spatial
distribution of VO measured by SPM can be successfully predicted and
interpreted by first-principles calculations if the coupling between VO 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 VO-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$\times$10$\,\text{nm}^{2}$) surface morphologies resembling the SPM
measurements. This analysis revealed new physical properties of the polarons
on TiO2(110), where the formation of inhomogeneously distributed VO is linked
to an increased formation of surface polarons and, therefore, to the density
of active sites.

##  2 Results and Discussion

###  2.1 Defect Distribution via DFT, Experiment, and Machine Learning

Figure 1: Oxygen vacancy distribution on rutile TiO2(110) obtained by various
methods. (a) Schematic representation of the most favorable VO distribution in
non-polaronic DFT calculations as obtained from a 6$\times$4
($\sim$1.8$\times$2.6$\,\text{nm}^{2}$) supercell. The schematic depiction is
generated by showing the Obr bridging atoms as black regions and Ti5c rows and
VO as white. The inset displays the structural model of rutile TiO2(110). The
distance maximizing VO distribution (6 sites in row, 3 sites in adjacent row)
and the 6$\times$4 supercell are indicated. (b) Unoccupied-states, constant-
current STM image of a clean, reduced rutile TiO2(110) surface (imaging
parameters in the Figure) depicting Ti5c rows and VOs as bright, while Obr
rows are depicted as dark. More details on the contrast formation are given in
the Methods. Locally low and high VO-concentrations (c${}_{\rm 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${}_{\textrm{S0}}$ and
Pol${}_{\textrm{S1}}$, respectively) and VOs are modeled in a 54$\times$24
($\sim$ 16$\times$16$\,\text{nm}^{2}$) supercell. Orange and yellow markers
show the position of surface and subsurface polarons in the ML prediction.

Fig. 1 shows the surface structure of reduced rutile TiO2(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$\times$1 rutile surface consists of
alternating rows of under-coordinated (two-fold) oxygen atoms (the bridging
oxygen atoms, Obr) and five-fold coordinated titanium atoms (Ti5c) running
along the [001] direction 38, 39. Oxygen vacancies form easily on the Obr
sites upon sputtering and annealing, up to a critical concentration of
c${}_{\rm V_{O}}$$\simeq 17$% 34. At stronger reducing conditions, the surface
undergoes a structural reconstruction40, 41, 42, 43, 44, 45. Every VO releases
two excess electrons that form polaronic states, localizing preferably on
subsurface Ti sites 27, 46, 35. Thus, the VO 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 VO-VO distances represents the most favorable vacancy
distribution. For the critical concentration of c${}_{\rm V_{O}}$$=17$%, this
corresponds to a homogeneous configuration with a VO-VO 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 VO 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 TiO2 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 27. While 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 28.

The experimental measurements do not support such homogeneous VO pattern. Fig.
1b shows a typical image as obtained from low-temperature STM measurements on
a TiO2(110) surface after sputtering and annealing treatment to form a high
content of oxygen vacancies (c${}_{\rm V_{O}}$$\simeq 14$%, close to the
critical value of $17$%). At this temperature, the oxygen vacancies (imaged as
bright spots along the dark [001] Obr 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 VO surface
structure. Simply adding polarons on a rigid VO pattern (effectively
decoupling VO 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 VO and polaron positions to
find the configuration that minimizes the total energy of the system. The
resulting VO 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 VO-VO interactions. Moreover, it allows us to
consider large surface areas of about 250$\,\text{nm}^{2}$
(>15$\times$15$\,\text{nm}^{2}$), corresponding to 54$\times$24 supercells,
extending considerably the limits of standard DFT simulations.

###  2.2 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$\times$4 and 12$\times$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 48.
Finally, 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 TiO2(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 52, 27. Polarons were
localized at chosen surface Ti${}_{\textrm{S0}}$ (Pol${}_{\textrm{S0}}$) and
subsurface Ti${}_{\textrm{S1}}$ (Pol${}_{\textrm{S1}}$) sites via occupation
matrix control 53. We modeled 2367 symmetrically-inequivalent polaron-VO
configurations in a 6$\times$4 supercell (i.e., six and four times the [001]
and [1$\bar{1}$0] lattice vectors, respectively), and 2155 configurations in a
12$\times$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: 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.

Fig. 2 sketches the main features of our ML architecture predicting the
stability of different defect-polaron configurations A generic polaron-VO
distribution on the rutile TiO2 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 27. To 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):

$E_{\textrm{tot}}=\sum_{i}^{N_{\textrm{S1}}}E_{i}+\sum_{j}^{N_{\textrm{S0}}}E_{j}+\sum_{k}^{N_{\textrm{VO}}}E_{k}~{}~{}\textrm{.}$
(1)

Here, $E_{\textrm{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${}_{\textrm{S1}}$, Pol${}_{\textrm{S0}}$, and VO 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 27. The 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
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$\times$4
unit cell (see SM Fig. S1), we achieved a mean absolute error of $1.8$ and
$2.2$ meV/VO for the training and validation sets, respectively. By adding
training data from the 12$\times$2 unit cell (see SM Fig. S2), the mean
absolute error increased slightly ($2.9$ and $3.5$ meV/VO 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\bar{1}0$] directions. For a detailed description of the training
process see Methods Section.

Aiming for a comparison with the experimental measurements, we focus here on
the low-energy 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, VO on rutile are immobile, while polarons always show a
certain degree of mobility, hopping/tunneling a few lattice sites around the
equilibrium position 39, 34. The simulated annealing can be implemented as a
global optimization scheme 54. Candidate 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$\times$10.5$\,\text{nm}^{2}$ to collect statistics
and 16$\times$16$\,\text{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$\times$24
rutile TiO2(110) supercell. In the following, we quantitatively analyze our
results.

###  2.3 Formation of VO-Polaron Patterns and Their Mutual Interaction

Figure 3: Analysis of VO-polaron configurations in a 6$\times$4 TiO2(110)
cell. (a) Comparison of selected low-energy vacancy-polaron configurations as
obtained by different treatments of the polaron-VO 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 $\Delta E$ are shown using the
homogeneous VO distributions (from “No Polarons” and “Random Polarons”) as
references (note the large energy gain $\Delta E_{\mathrm{pol}}$ of -3.23 eV
between non-polaronic and polaronic solutions with homogeneous VO patterns).
Dashed lines connect identical VO configurations. New VO configurations found
in the ML search are displayed in red. The occurrence of Ti${}_{\textrm{S0}}$
polarons is highlighted in orange and purple for [1$\bar{1}$2]- and
[1$\bar{1}$0]-aligned oxygen vacancies, respectively. The most important VO-
polaron complexes are shown schematically in top view at the bottom of the
Figure. Only the most stable polaronic configuration per VO arrangement is
shown. (b,c) Top and side views of the polaronic isocharge surfaces of the
[1$\bar{1}$2]-aligned VO-Pol${}_{\textrm{S0}}$ complex (top), and of the
Pol${}_{\textrm{S1}}$ in the homogeneous VO-distribution (bottom).

The analysis of the low-energy configurations (see SM Fig. S4 for the energy
distribution of all possible VO configurations in the 6$\times$4 cell) is
summarized in Fig. 3. Fig. 3a shows the improvement of energies of the
TiO2(110) 6$\times$4 cell as obtained by treating VO-polaron coupling at three
different levels: (i) Suppressing polarons (“No Polarons”); (ii) Distributing
polarons in random or positions guided by physical intuition (“Random
Polarons” ); (iii) Full inclusion of polaron-VO 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${}_{\rm NoP}^{\rm Hom}$ 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 37, 34, 27. Here, the ground state configuration
is given by a homogeneous distribution of Pol${}_{\textrm{S1}}$ in the
homogeneous VO pattern (labeled as C${}_{\rm RandP}^{\rm Hom}$ in Fig. 3a).

Treating polaron-VO coupling at the ML level (third approach) results in novel
VO 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${}_{\textrm{S1}}$ showing a better
stability in the homogeneous VO background (labeled as C${}_{\rm ML}^{\rm
Hom}$, see also SM Fig. S5). Moreover, new polaron configurations explored by
the extensive ML search improve the stability of many other VO 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${}_{\rm ML}^{\rm 0}$ and
C${}_{\rm ML}^{\rm 2}$) resulted in energy values even lower than the
homogeneous distribution, revealing a new ground state for the system.
Moreover, novel VO-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${}_{\rm ML}^{\rm 1}$). The polarons play a
key role in stabilizing this VO-pattern and as further proof, we calculated
the energy of this new VO-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${}_{\rm ML}^{\rm Hom}$), we note the presence of
at least one polaron on a surface Ti${}_{\textrm{S0}}$ site (configurations
containing Pol${}_{\textrm{S0}}$ are orange highlighted in Fig. 3a). Fig. 3b
and c compare the spatial distribution of the surface Pol${}_{\textrm{S0}}$
and subsurface Pol${}_{\textrm{S1}}$. The formation of the surface polaron is
particularly stable when occurring in the central Ti${}_{\textrm{S0}}$ site
between two oxygen vacancies aligned on the [1$\bar{1}$2] direction (see top
view in Fig. 3b). This [1$\bar{1}$2]-aligned VO-polaron complex represents
indeed the ground state configuration obtained by our ML search (e.g., it is
present in C${}_{\rm ML}^{\rm 0,1,2}$). Another remarkably stable complex is
given by two vacancies aligned along the [$1\bar{1}0$] direction and one
Pol${}_{\textrm{S0}}$ in their vicinity (see SM Fig. S5). For instance, this
complex appears in the configurations highlighted in purple in Fig. 3a
(C${}_{\rm ML}^{\rm 10,11}$). The [1$\bar{1}$2] and [$1\bar{1}0$] 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 VO-
regions. In contrast, DFT predictions, which neglect polaron-VO interaction or
randomly distribute polarons, favor homogeneous configurations.

###  2.4 Comparison of a Large Scale Model and the Experimental Surface

Figure 4: Autocorrelation functions of the VO 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 VO autocorrelation functions
along a single (a) and adjacent (b) [001]-aligned Obr 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$\times$16 supercells. ML Pol and ML No Pol are started from identically
seeded VO patterns, with 2$N_{\mathrm{V_{O}}}$ polarons and no polarons,
respectively. Autocorrelation functions are rescaled to account for c${}_{\rm
V_{O}}$ differences in experiment (14.2%) and simulation (16.7%).

Fig. 4 shows our results as obtained by ML-driven annealing simulations on
large-area 24$\times$16 cells (corresponding to 7$\times$10$\,\text{nm}^{2}$),
which enables a direct comparison with the experiment. Visual inspection (Fig.
1b vs. c) already indicates that our ML treatment provides VO 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 VO distributions extracted from Fig. 1b (for details see Supplementary
Fig. S7). The simulated annealing procedure starts from random VO-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 VO pattern but takes into account
annealing-induced disorder effects.

The ACFs are shown in Fig. 4, where projections of VO 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 VO-VO 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-VO
interactions (see Fig. 4a and b red and blue data, respectively). By
considering only the VO-VO 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-VO 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 VO arrangement.

As a result, the rutile TiO2 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$\bar{1}$2]- and
[$1\bar{1}0$]-aligned VO-polaron complexes contribute to this alternation of
locally less and more reduced areas at a given c${}_{\rm V_{O}}$. To further
corroborate this result, we performed additional DFT calculations modeling
this strong inhomogeneity (see SM Fig. S5 with configurations C${}_{\rm
ML}^{\rm 26-29}$). We also note that this analysis reconciles the DFT
predictions on the critical concentration at which the (1$\times$2) surface
reconstruction occurs for the surface phase transition, which was calculated
as $\sim 20\%$, in apparent disagreement with the experiments reporting an
average concentration of $17$% 34.

##  3 Conclusions

In summary, we directly elucidated the impact of polarons on the structure of
oxide surfaces, using an example of the prototypical rutile TiO2(110) surface.
Specifically, we designed a computational machinery to predict the
distribution of polarons and oxygen vacancies on rutile TiO2(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 VO distribution as detected by the experiments. Our
analysis clarifies the peculiar inhomogeneous distribution of VO on rutile
TiO2(110). Most importantly, the system shows a tendency towards the formation
of high-density VO patterns, alternated with low-density VO regions. While
larger defect-free areas are typically attributed to subsurface Ar impurities
58, the here observed fluctuation of the local c${}_{\rm V_{O}}$ can partially
be attributed to the interaction of polarons and VOs.

These results suggest that surface reactivity could be optimized by tuning the
annealing procedure to facilitate the formation of energetically more
favorable, high-density VO patterns, which promote surface localized charges
and their interaction with adsorbates 26. To 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${}_{\rm V_{O}}$59. SPM measurements in the presence of CO adsorbates
confirm the formation of the [1$\bar{1}$2]-aligned VO-polaron complex 26,
while STM measurements probing the filled states on the clean surface do show
some disparity in comparison to simulated STM37. The 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
SrTiO3(001) surface 61 exhibiting Sr adatom/vacancy and often doped by Nb
atoms 29. Moreover, 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.

##  4 Methods

###  4.1 DFT Modeling

We performed DFT+U calculations using VASP 49, 50, 51 on the rutile TiO2(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 52, 27. The
sampling of the reciprocal space included the $\Gamma$-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$\times$4 and 12$\times$2\. To partially account for the role of
thermal effects in the stabilization of the VO 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 measurements62, 63. This strain of $+0.5$% is well below the
crossover point of $+3$%, where surface polaron formation is favored over
subsurface polaron formation60.

Within the supercells, we removed 4 surface bridging oxygen atoms (in random
positions) from every slab, obtaining a c${}_{\rm V_{O}}$ of approximately
$17$%. To assess non-polaronic solutions we performed spin unpolarized DFT,
constraining the excess electrons in spatially delocalized 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${}_{\textrm{S1}}$ and Pol${}_{\textrm{S0}}$ sites
27. Finally, we performed an unconstrained relaxation starting from the
structures and wave functions determined in the previous step.

###  4.2 ML Model Training and Optimization of Defect Configurations

The machine learning model is implemented in the framework of JAX 64. Here, we
describe the model optimization based on the study of configurations in the
6$\times$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 VO-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.

###  4.3 Experimental Setup

SPM was performed using STM in an ultrahigh vacuum (UHV) chamber with a base
pressure below 2$\times 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. Stiff qPlus sensors 67 ($k=1800$ N$\cdot$m-1,
$Q$=5000–30000, $f_{0}$ $\in$ [25-45] kHz) with a a sharp W tip 68 were used
to collect the tunneling current ($I_{\mathrm{t}}$) and the frequency shift
($\Delta f$) signals; deflection detection was achieved using a cryogenic
preamplifier in vacuum 69. W 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 TiO2(110) surface. Tip sharpness was indicated by the
low frequency shifts ($\Delta f$ $\in$ $[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$\times$10-10 mbar. Surfaces were cleaned by cycles of
sputtering and UHV annealing that consequently reduced the samples and
introduced VOs 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
$\mu$A$\cdot$ 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$\times 10^{-7}$ mbar of O2 shower for 10 min. When the
reduction level was too high, the rutile TiO2 samples were re-oxidized $ex-
situ$ at 800 ∘C in O2 flow and reintroduced to UHV for cleaning.

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_{\mathrm{S}}$ = +0.9 V, grounded tip, tunneling current set-point
$I_{\mathrm{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 TiO2(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 VO states, while the
O states constitute the valence band 70. Therefore, the highest probability of
electron tunneling from the tip to the surface is above the Ti5c rows and VOs
– they appear bright under these STM conditions. On the other hand, the
tunneling is less likely above the Obr rows and they appear dark even though
they geometrically protrude highest from the surface. Note that in Fig. 1b
individual Ti5c atoms can be recognized as spheres forming a row along the
$[001]$ direction, while VOs are recognized as isolated, bright spheres.

##  5 Data availability

The data presented in this article is available with the accompanying code
upon publication or from the corresponding author upon request.

##  6 Code availability

A minimal example of the code used to produce the presented results is
available in a Github repository upon publication.

##  7 Author Contributions

VB implemented the ML model and performed calculations together with MP and
MR. VB wrote the first draft under the supervision of CF and MR. CF and MR
conceptualized and supervised the work. IS and MS performed experiments under
the supervision of UD. All authors contributed to reviewing and editing the
final draft.

##  8 Competing Interests

The authors declare no competing interests.

This work was supported by the Austrian Science Fund (FWF) project SFB-F81
project TACO. The computational results have been achieved using the Vienna
Scientific Cluster (VSC). CF acknowledges the NextGenerationEU-Piano Nazionale
Resistenza e Resilienza (PNRR) CN-HPC grant no. (CUP) J33C22001170001, SPOKE 6
- Multiscale Modelling & Engineering. KP acknowledges the NRDIO-Hungary grant
no. FK124100 and a Bolyai Fellowship of the Hungarian Academy of Sciences.

## References

  * Rousseau et al. 2020 Rousseau, R.; Glezakou, V.-A.; Selloni, A. Theoretical insights into the surface physics and chemistry of redox-active oxides. _Nature Reviews Materials_ 2020, _5_ , 460–475
  * Franceschi and Diebold 2023 Franceschi, G.; Diebold, U. _Encyclopedia of Materials: Electronics_ ; 2023; pp 501–511
  * Jupille and Thornton 2015 Jupille, J.; Thornton, G. _Defects at Oxide Surfaces_ ; Springer Series in Surface Sciences; 2015; Vol. 58; pp 327–349
  * Strand and Shluger 2023 Strand, J.; Shluger, A. L. On the Structure of Oxygen Deficient Amorphous Oxide Films. _Advanced Science_ 2023, 2306243
  * Franchini et al. 2021 Franchini, C.; Reticcioli, M.; Setvin, M.; Diebold, U. Polarons in materials. _Nature Reviews Materials_ 2021, _6_ , 560–586
  * Emin 2013 Emin, D. _Polarons_ ; Cambridge University Press, 2013
  * Alexandrov and Devreese 2010 Alexandrov, A. S.; Devreese, J. T. _Advances in Polaron Physics_ ; Springer Series in Solid-State Sciences; 2010; Vol. 159; p 171
  * Stoneham et al. 2007 Stoneham, A. M.; Gavartin, J.; Shluger, A. L.; Kimmel, A. V.; Mũoz Ramo, D.; Rønnow, H. M.; Aeppli, G.; Renner, C. Trapping, self-trapping and the polaron family. _Journal of Physics Condensed Matter_ 2007, _19_ , 255208
  * Pastor et al. 2022 Pastor, E.; Sachs, M.; Selim, S.; Durrant, J. R.; Bakulin, A. A.; Walsh, A. Electronic defects in metal oxide photocatalysts. _Nature Reviews Materials_ 2022, _7_ , 503–521
  * Kick et al. 2020 Kick, M.; Grosu, C.; Schuderer, M.; Scheurer, C.; Oberhofer, H. Mobile Small Polarons Qualitatively Explain Conductivity in Lithium Titanium Oxide Battery Electrodes. _The Journal of Physical Chemistry Letters_ 2020, _11_ , 2535–2540
  * Chen et al. 2023 Chen, M.; Grieder, A. C.; Smart, T. J.; Mayford, K.; McNair, S.; Pinongcos, A.; Eisenberg, S.; Bridges, F.; Li, Y.; Ping, Y. The impacts of dopants on the small polaron mobility and conductivity in hematite – the role of disorder. _Nanoscale_ 2023, _15_ , 1619–1628
  * Smart and Ping 2017 Smart, T. J.; Ping, Y. Effect of defects on the small polaron formation and transport properties of hematite from first-principles calculations. _Journal of Physics Condensed Matter_ 2017, _29_ , 394006
  * Cheng et al. 2022 Cheng, C.; Zhu, Y.; Fang, W.-H.; Long, R.; Prezhdo, O. V. CO Adsorbate Promotes Polaron Photoactivity on the Reduced Rutile TiO2(110) Surface. _JACS Au_ 2022, _2_ , 234–245
  * Cheng et al. 2022 Cheng, C.; Zhu, Y.; Zhou, Z.; Long, R.; Fang, W.-H. Photoinduced small electron polarons generation and recombination in hematite. _npj Computational Materials_ 2022, _8_ , 1–8
  * Sokolović et al. 2020 Sokolović, I.; Reticcioli, M.; Čalkovský, M.; Wagner, M.; Schmid, M.; Franchini, C.; Diebold, U.; Setvín, M. Resolving the adsorption of molecular O2 on the rutile TiO2(110) surface by noncontact atomic force microscopy. _Proceedings of the National Academy of Sciences of the United States of America_ 2020, _117_ , 14827–14837
  * Tanner et al. 2021 Tanner, A. J.; Wen, B.; Ontaneda, J.; Zhang, Y.; Grau-Crespo, R.; Fielding, H. H.; Selloni, A.; Thornton, G. Polaron-Adsorbate Coupling at the TiO2(110)-Carboxylate Interface. _The Journal of Physical Chemistry Letters_ 2021, _12_ , 3571–3576
  * Yim et al. 2018 Yim, C. M.; Chen, J.; Zhang, Y.; Shaw, B.-J.; Pang, C. L.; Grinter, D. C.; Bluhm, H.; Salmeron, M.; Muryn, C. A.; Michaelides, A.; Thornton, G. Visualization of Water-Induced Surface Segregation of Polarons on Rutile TiO2(110). _The Journal of Physical Chemistry Letters_ 2018, _9_ , 4865–4871
  * Cheng et al. 2023 Cheng, C.; Zhou, Z.; Long, R. Time-Domain View of Polaron Dynamics in Metal Oxide Photocatalysts. _The Journal of Physical Chemistry Letters_ 2023, _14_ , 10988–10998
  * Ren et al. 2023 Ren, Z.; Shi, Z.; Feng, H.; Xu, Z.; Hao, W. Recent Progresses of Polarons: Fundamentals and Roles in Photocatalysis and Photoelectrocatalysis. _Advanced Science_ 2023, 2305139
  * Dohnálek et al. 2010 Dohnálek, Z.; Lyubinetsky, I.; Rousseau, R. Thermally-driven processes on rutile TiO2(110)-(1×1): A direct view at the atomic scale. _Progress in Surface Science_ 2010, _85_ , 161–205
  * Tanner and Thornton 2022 Tanner, A. J.; Thornton, G. TiO2 Polarons in the Time Domain: Implications for Photocatalysis. _The Journal of Physical Chemistry Letters_ 2022, _13_ , 559–566
  * Sombut et al. 2022 Sombut, P.; Puntscher, L.; Atzmueller, M.; Jakub, Z.; Reticcioli, M.; Meier, M.; Parkinson, G. S.; Franchini, C. Role of Polarons in Single-Atom Catalysts: Case Study of Me1 [Au1, Pt1, and Rh1] on TiO2(110). _Topics in Catalysis_ 2022, _2_ , 1–16
  * Geiger and López 2022 Geiger, J.; López, N. Coupling Metal and Support Redox Terms in Single-Atom Catalysts. _The Journal of Physical Chemistry C_ 2022, _126_ , 13698–13704
  * Geiger et al. 2022 Geiger, J.; Sabadell-Rendón, A.; Daelman, N.; López, N. Data-driven models for ground and excited states for Single Atoms on Ceria. _npj Computational Materials_ 2022, _8_ , 1–8
  * Cao et al. 2017 Cao, Y.; Yu, M.; Qi, S.; Huang, S.; Wang, T.; Xu, M.; Hu, S.; Yan, S. Scenarios of polaron-involved molecular adsorption on reduced TiO2(110) surfaces. _Scientific Reports_ 2017, _7_ , 6148
  * Reticcioli et al. 2019 Reticcioli, M.; Sokolović, I.; Schmid, M.; Diebold, U.; Setvin, M.; Franchini, C. Interplay between Adsorbates and Polarons: CO on Rutile TiO ${}_{\textrm{2}}$(110). _Physical Review Letters_ 2019, _122_ , 016805
  * Birschitzky et al. 2022 Birschitzky, V. C.; Ellinger, F.; Diebold, U.; Reticcioli, M.; Franchini, C. Machine learning for exploring small polaron configurational space. _npj Computational Materials_ 2022, _8_ , 1–9
  * Zhang et al. 2019 Zhang, D.; Han, Z. K.; Murgida, G. E.; Ganduglia-Pirovano, M. V.; Gao, Y. Oxygen-Vacancy Dynamics and Entanglement with Polaron Hopping at the Reduced CeO2 (111) Surface. _Physical Review Letters_ 2019, _122_ , 096101
  * Ellinger et al. 2023 Ellinger, F.; Shafiq, M.; Ahmad, I.; Reticcioli, M.; Franchini, C. Small Polaron Formation on the Nb-doped SrTiO3(001) Surface. _Physical Review Materials_ 2023, _7_ , 064602
  * Österbacka et al. 2022 Österbacka, N.; Ambrosio, F.; Wiktor, J. Charge Localization in Defective BiVO${}_{\textrm{4}}$. _The Journal of Physical Chemistry C_ 2022, _126_ , 2960–2970
  * Sun et al. 2017 Sun, L.; Huang, X.; Wang, L.; Janotti, A. Disentangling the role of small polarons and oxygen vacancies in CeO2. _Physical Review B_ 2017, _95_ , 245101
  * Reticcioli et al. 2019 Reticcioli, M.; Diebold, U.; Kresse, G.; Franchini, C. _Handbook of Materials Modeling_ ; Springer International Publishing, 2019; pp 1–39
  * Pham and Deskins 2020 Pham, T. D.; Deskins, N. A. Efficient Method for Modeling Polarons Using Electronic Structure Methods. _Journal of Chemical Theory and Computation_ 2020, _16_ , 5264–5278
  * Reticcioli et al. 2017 Reticcioli, M.; Setvin, M.; Hao, X.; Flauger, P.; Kresse, G.; Schmid, M.; Diebold, U.; Franchini, C. Polaron-driven surface reconstructions. _Physical Review X_ 2017, _7_ , 031053
  * Kowalski et al. 2010 Kowalski, P. M.; Camellone, M. F.; Nair, N. N.; Meyer, B.; Marx, D. Charge localization dynamics induced by oxygen vacancies on the TiO2(110) surface. _Physical Review Letters_ 2010, _105_ , 146405
  * Han et al. 2018 Han, Z. K.; Yang, Y. Z.; Zhu, B.; Ganduglia-Pirovano, M. V.; Gao, Y. Unraveling the oxygen vacancy structures at the reduced CeO2(111) surface. _Physical Review Materials_ 2018, _2_ , 035802
  * Reticcioli et al. 2018 Reticcioli, M.; Setvin, M.; Schmid, M.; Diebold, U.; Franchini, C. Formation and dynamics of small polarons on the rutile TiO2(110) surface. _Physical Review B_ 2018, _98_ , 045306
  * Diebold 2003 Diebold, U. The surface science of titanium dioxide. _Surface Science Reports_ 2003, _48_ , 53–229
  * Setvin et al. 2014 Setvin, M.; Franchini, C.; Hao, X.; Schmid, M.; Janotti, A.; Kaltak, M.; Van De Walle, C. G.; Kresse, G.; Diebold, U. Direct view at excess electrons in TiO2 rutile and anatase. _Physical Review Letters_ 2014, _113_ , 086402
  * Onishi and Iwasawa 1994 Onishi, H.; Iwasawa, Y. Reconstruction of TiO2(110) surface: STM study with atomic-scale resolution. _Surface Science_ 1994, _313_ , L783–L789
  * Li et al. 1999 Li, M.; Hebenstreit, W.; Gross, L.; Diebold, U.; Henderson, M. A.; Jennison, D. R.; Schultz, P. A.; Sears, M. P. Oxygen-induced restructuring of the TiO2(110) surface: a comprehensive study. _Surface Science_ 1999, _437_ , 173–190
  * Li et al. 2000 Li, M.; Hebenstreit, W.; Diebold, U. Morphology change of oxygen-restructured surfaces by UHV annealing: Formation of a low-temperature structure. _Physical Review B_ 2000, _61_ , 4926–4933
  * McCarty and Bartelt 2003 McCarty, K. F.; Bartelt, N. C. The $1\times 1$/$1\times 2$ phase transition of the TiO2(110) surface-variation of transition temperature with crystal composition. _Surface Science_ 2003, _527_ , L203–L212
  * Wang et al. 2014 Wang, Q.; Oganov, A. R.; Zhu, Q.; Zhou, X. F. New reconstructions of the (110) surface of rutile TiO2 predicted by an evolutionary method. _Physical Review Letters_ 2014, _113_ , 266101
  * Mochizuki et al. 2016 Mochizuki, I.; Ariga, H.; Fukaya, Y.; Wada, K.; Maekawa, M.; Kawasuso, A.; Shidara, T.; Asakura, K.; Hyodo, T. Structure determination of the rutile-TiO2(110)-(1$\times$2) surface using total-reflection high-energy positron diffraction (TRHEPD). _Physical Chemistry Chemical Physics_ 2016, _18_ , 7085–7092
  * Deskins et al. 2009 Deskins, N. A.; Rousseau, R.; Dupuis, M. Localized electronic states from surface hydroxyls and polarons in TiO2(110). _Journal of Physical Chemistry C_ 2009, _113_ , 14583–14586
  * LeCun et al. 2015 LeCun, Y.; Bengio, Y.; Hinton, G. Deep learning. _Nature_ 2015, _521_ , 436–444
  * Behler 2021 Behler, J. Four Generations of High-Dimensional Neural Network Potentials. _Chemical Reviews_ 2021, _121_ , 10037–10072
  * Kresse and Furthmüller 1996 Kresse, G.; Furthmüller, J. Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set. _Computational Materials Science_ 1996, _6_ , 15–50
  * Kresse and Joubert 1999 Kresse, G.; Joubert, D. From ultrasoft pseudopotentials to the projector augmented-wave method. _Physical Review B_ 1999, _59_ , 1758–1775
  * Blöchl 1994 Blöchl, P. E. Projector augmented-wave method. _Physical Review B_ 1994, _50_ , 17953–17979
  * Dudarev and Botton 1998 Dudarev, S.; Botton, G. Electron-energy-loss spectra and the structural stability of nickel oxide: An LSDA+U study. _Physical Review B_ 1998, _57_ , 1505–1509
  * Allen and Watson 2014 Allen, J. P.; Watson, G. W. Occupation matrix control of d- and f-electron localisations using DFT + U. _Physical Chemistry Chemical Physics_ 2014, _16_ , 21016–21031
  * Kirkpatrick et al. 1983 Kirkpatrick, S.; Gelatt, C. D.; Vecchi, M. P. Optimization by Simulated Annealing. _Science_ 1983, _220_ , 671–680
  * Metropolis et al. 1953 Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M. N.; Teller, A. H.; Teller, E. Equation of State Calculations by Fast Computing Machines. _The Journal of Chemical Physics_ 1953, _21_ , 1087–1092
  * Ji et al. 2022 Ji, W.; Wang, N.; Li, Q.; Zhu, H.; Lin, K.; Deng, J.; Chen, J.; Zhang, H.; Xing, X. Oxygen vacancy distributions and electron localization in a CeO2(100) nanocube. _Inorganic Chemistry Frontiers_ 2022, _9_ , 275–283
  * Setvin et al. 2015 Setvin, M. et al. A Multitechnique Study of CO Adsorption on the TiO2 Anatase (101) Surface. _Journal of Physical Chemistry C_ 2015, _119_ , 21044–21052
  * Potapenko et al. 2014 Potapenko, D. V.; Li, Z.; Kysar, J. W.; Osgood, R. M. Nanoscale Strain Engineering on the Surface of a Bulk TiO2 Crystal. _Nano Letters_ 2014, _14_ , 6185–6189
  * Krüger et al. 2008 Krüger, P.; Bourgeois, S.; Domenichini, B.; Magnan, H.; Chandesris, D.; Le Fèvre, P.; Flank, A. M.; Jupille, J.; Floreano, L.; Cossaro, A.; Verdini, A.; Morgante, A. Defect states at the TiO2(110) surface probed by resonant photoelectron diffraction. _Physical Review Letters_ 2008, _100_ , 055501
  * Reticcioli et al. 2022 Reticcioli, M.; Diebold, U.; Franchini, C. Modeling polarons in density functional theory: lessons learned from TiO2. _Journal of Physics: Condensed Matter_ 2022, _34_ , 204006
  * Sokolović et al. 2019 Sokolović, I.; Schmid, M.; Diebold, U.; Setvín, M. Incipient ferroelectricity: A route towards bulk-terminated SrTiO3. _Phys. Rev. Mater._ 2019, _3_ , 034407
  * Kirby 1967 Kirby, R. K. Thermal Expansion of Rutile from 100 to 700 °K. _Journal of Research of the National Bureau of Standards. Section A, Physics and Chemistry_ 1967, _71A_ , 363–369
  * Hummer et al. 2007 Hummer, D. R.; Heaney, P. J.; Post, J. E. Thermal expansion of anatase and rutile between 300 and 575 K using synchrotron powder X-ray diffraction. _Powder Diffraction_ 2007, _22_ , 352–357
  * Bradbury et al. 2018 Bradbury, J.; Frostig, R.; Hawkins, P.; Johnson, M. J.; Leary, C.; Maclaurin, D.; Necula, G.; Paszke, A.; VanderPlas, J.; Wanderman-Milne, S.; Zhang, Q. JAX: composable transformations of Python+NumPy programs. 2018; http://github.com/google/jax
  * Zhang et al. 2007 Zhang, Z.; Ge, Q.; Li, S.-C.; Kay, B. D.; White, J. M.; Dohnálek, Z. Imaging Intrinsic Diffusion of Bridge-Bonded Oxygen Vacancies on TiO2( 110 ). _Physical Review Letters_ 2007, _99_ , 126105
  * Schmid et al. 2019 Schmid, M.; Setvín, M.; Diebold, U. Device for suspending a load in a vibration-insulated manner. 2019
  * Giessibl 2019 Giessibl, F. J. The qPlus sensor, a powerful core for the atomic force microscope. _Rev. Sci. Instr._ 2019, _90_ , 011101
  * Setvín et al. 2012 Setvín, M.; Javorskỳ, J.; Turčinková, D.; Matolínová, I.; Sobotík, P.; Kocán, P.; Ošt’ádal, I. Ultrasharp tungsten tips—characterization and nondestructive cleaning. _Ultramicroscopy_ 2012, _113_ , 152–157
  * Huber and Giessibl 2017 Huber, F.; Giessibl, F. J. Low noise current preamplifier for qPlus sensor deflection signal detection in atomic force microscopy at room and low temperatures. _Rev. Sci. Instrum._ 2017, _88_ , 073702
  * Diebold et al. 1996 Diebold, U.; Anderson, J. F.; Ng, K. O.; Vanderbilt, D. Evidence for the Tunneling Site on Transition-Metal Oxides: TiO2(110). _Phys. Rev. Lett._ 1996, _77_ , 1322