Modeling polarons in density functional theory: lessons learned from TiO₂

Excess charge in polarizable materials can couple to the lattice phonons and form polaron quasiparticles [1]. Since their formulation in the early twentieth century [2–5], polarons have extensively been studied in physics, materials science and chemistry [6–10], and continue nowadays to be of large interest due also to their impact on a wide range of modern technologies [11–16]. Titanium dioxide is a prototypical example for the study of polarons in materials [16–21]. In its rutile form, TiO₂ is likely to host so-called small electron polarons, i.e., electrons strongly localized on one or few lattice sites, and associated to evident distortions of the local structure [22, 23]. The long series of diverse experiments performed on this compound [1, 21] has been accompanied by analysis at the density-functional theory (DFT) level, which is a broadly used technique for investigations of polaronic properties [24–39]. Modelization of the localized electrons requires a correction of the electronic correlation of the Ti d orbitals, that can be efficiently achieved by adopting the DFT + U technique: an additional energy term is introduced in the Hamiltonian, and tuned by a parameter U [40]. The value of U is typically determined by comparing the electronic properties of the simulated system with the experimental measurements, or by relying on different ab initio approaches [41–46].

Figure 1 collects general aspects of polaron formation on the rutile TiO₂(110) surface as obtained by DFT + U calculations. Oxygen vacancies (V_O) that easily develop on the surface of rutile samples introduce excess electrons into the system. Assuming a rigid band model, one would expect this excess charge to fill the bottom of the conduction band, with a transition from a semiconducting density of states (DOS), typical of the pristine surface (figure 1(b)), to a metallic phase (figure 1(c)): we refer to this system as delocalized solution, since the electronic states are spatially dispersed throughout the material (figure 1(e)). Formation of a polaronic state with the excess charge localized at single Ti atoms and associated to a sharp state in the energy gap (figures 1(a) and (d)) is predicted by DFT + U as a more stable solution. The charge localization can occur at different lattice sites, leading to polarons characterized by different properties and stability: usually, the most suitable hosting sites are the subsurface Ti atoms lying below the five-fold coordinated atoms (${\mathrm{T}\mathrm{i}}_{S1}^{A}$, see labeling in figure 1) [47–50]. However, environmental conditions and/or computational details can alter the polaron stability, and favor the charge localization elsewhere, especially on the surface ${\mathrm{T}\mathrm{i}}_{S0}^{A}$ atoms [26, 27, 43, 44, 51–55].

Zoom In Zoom Out Reset image size

Here, we systematically analyze the polaronic properties of rutile TiO₂(110), and the dependence on physical quantities, such as lattice strain, and the correction U of the electronic correlation in DFT + U calculations. In particular, we describe the change in the charge localization, local structural distortions and electronic properties of polarons as determined by different values of the U parameter. Moreover, we show that the application of strain along the [001] direction can reverse the stability of ${\mathrm{T}\mathrm{i}}_{S1}^{A}$ and ${\mathrm{T}\mathrm{i}}_{S0}^{A}$ polarons, in agreement with experimental observations, due to their different orbital symmetries and different interaction with the surrounding local structure. We emphasize that a precise knowledge of the effects of the parameters used in DFT calculations is key to obtain reliable predictions, in agreement with experiments: we conclude the manuscript by showing the effects of the altered relative stability of ${\mathrm{T}\mathrm{i}}_{S1}^{A}$ and ${\mathrm{T}\mathrm{i}}_{S0}^{A}$ polarons on a well studied application, the adsorption of CO molecules on the rutile surface [56].

This study was conducted in the framework of density functional theory (DFT), by using the Vienna ab initio simulation package (VASP) [57–59]. We adopted the generalized gradient approximation (GGA) with the Perdew, Burke, and Ernzerhof parametrization [60], optimized to account for Van-der-Waals dispersion interactions as originally proposed by Dion et al. (optPBE) [61, 62] in order to properly model the interaction with adsorbates. A correction for the electronic correlation was included by using an on-site effective U [40] of U^cRPA = 3.9 eV on the d orbitals of the Ti atoms [63] (optPBE + U calculations), previously determined by constrained random-phase approximation (cRPA) calculations in bulk rutile [35]. The rutile surfaces were modeled with asymmetric slabs containing five Ti layers in large two-dimensional 9 × 2 and 6 × 2 unit cells, using low-temperature experimental lattice parameters [21], with the inclusion of a vacuum region of 15 Å. In our reference frame, the x, y and z Cartesian axis correspond to the [001], [1$\bar{1}$0] and [110] crystallographic directions, respectively. All atomic sites except atoms on the bottom two layers were relaxed using standard convergence criteria with a plane-wave energy cutoff of 400 eV, and using the Γ point for sampling the reciprocal space.

Oxygen vacancies (V_O) were included by removing one two-fold coordinated surface O atom from the 9 × 2 slabs, and one or two O atoms from the 6 × 2 slabs, corresponding to concentrations of ${c}_{{V}_{\text{O}}}=$ 5.6%, 8.3%, and 16.7%, respectively. Every V_O supplies two excess electrons, eligible to form two polarons [35]. We used VESTA [64] to show the spatial extension of the excess charge. The polaron formation energy E_POL (${\mathrm{T}\mathrm{i}}_{S1}^{A}$) for sub-surface ${\mathrm{T}\mathrm{i}}_{S1}^{A}$ polarons distributed in equivalent localization sites in the slab is defined as the energy difference between the free energy of the system hosting polarons (E_loc) as compared to the energy of the system with the excess charge delocalized (E_del), divided by the number n of polarons in the slab: E_POL = (E_loc − E_del)/n. Single surface ${\mathrm{T}\mathrm{i}}_{S0}^{A}$ polarons were modeled by replacing one of the sub-surface ${\mathrm{T}\mathrm{i}}_{S1}^{A}$ polarons in the slab: the polaron formation energy E_POL (${\mathrm{T}\mathrm{i}}_{S0}^{A}$) was calculated by further subtracting the E_POL (${\mathrm{T}\mathrm{i}}_{S1}^{A}$) of the n_S1 sub-surface polarons, E_POL (${\mathrm{T}\mathrm{i}}_{S0}^{A}$) = E_loc − E_del − n_S1 E_POL (${\mathrm{T}\mathrm{i}}_{S1}^{A}$). The results for E_POL (${\mathrm{T}\mathrm{i}}_{S0}^{A}$) are comparable to the values obtained by modeling one single ${\mathrm{T}\mathrm{i}}_{S0}^{A}$ polaron and removing all remaining excess electron from the system, with the addition of a compensating background charge (results obtained in this setup are not shown here) [65]. The convergence of E_POL for the various polarons with respect to the modeling setup (number of atomic layers and lateral size of the slab) was discussed in previous works [65]. The CO adsorption energy (E_ads) is calculated by subtracting from the energy ${E}_{i}^{\text{CO}}$ of the system with the CO molecule adsorbing in a specific configuration i, the energy E_S1 of the clean slab hosting two ${\mathrm{T}\mathrm{i}}_{S1}^{A}$ polarons in their optimal configuration and the energy ${E}_{\text{gas}}^{\text{CO}}$ of the CO molecule in the gas phase [54]: ${E}_{\text{ads}}\ ={E}_{i}^{\text{CO}}-{E}_{S1}-{E}_{\text{gas}}^{\text{CO}}$.

Oxygen vacancies on the TiO₂(110) surface cause the formation of polarons, as described in the Introduction and shown in figure 1. Typically, sub-surface ${\mathrm{T}\mathrm{i}}_{S1}^{A}$ ions are the most favorable sites for polaron localization, and they have been described in detail in the literature [16–21, 27–29, 35, 47–50, 52–55]. These ${\mathrm{T}\mathrm{i}}_{S1}^{A}$ polarons show a ${d}_{{z}^{2}}$–${d}_{{x}^{2}-{y}^{2}}$ orbital symmetry, with the polaronic charge extending on the plane of the equatorial oxygen atoms, and hybridizing with surrounding Ti atoms on the same plane.

We note that on the surface ${\mathrm{T}\mathrm{i}}_{S0}^{A}$ sites, polaron formation occurs instead with different characteristics. Figure 2 and table 1 report the different solutions as obtained for a polaron on the next nearest neighbor ${\mathrm{N}\mathrm{N}\mathrm{N}-\mathrm{T}\mathrm{i}}_{S0}^{A}$ site from the V_O. The d_xz –d_yz symmetry in figures 2(a) and (b) is the most stable solution (i.e., a d_xz orbital, rotated from the [001] towards the direction of the bonds with two equatorial bonds), but the polaron can be stabilized also with the ${d}_{{x}^{2}-{y}^{2}}$ (figures 2(c) and (d)) and d_yz (figures 2(e) and (f)) symmetries, with similar energy stability (ΔE = +25 and +40 meV with respect to the d_xz –d_yz solution, respectively). These different solutions can be spontaneously obtained by simply adopting different starting conditions in our DFT simulations, driving the system towards different minima of the energy configuration space. Polarons lying far from the vacancy, seem to adopt preferentially a d_xz –d_yz symmetry, while the nearest neighbor ${\mathrm{N}\mathrm{N}-\mathrm{T}\mathrm{i}}_{S0}^{A}$ site to the vacancy allows for the stabilization of a ${d}_{{x}^{2}-{y}^{2}}$ symmetry, but these localization sites are in general less favorable than the ${\mathrm{N}\mathrm{N}\mathrm{N}-\mathrm{T}\mathrm{i}}_{S0}^{A}$ site [50, 65].

Zoom In Zoom Out Reset image size

Table 1. The ${\mathrm{N}\mathrm{N}\mathrm{N}-\mathrm{T}\mathrm{i}}_{S0}^{A}$ polaron. Orbital characters and relative energy stability (ΔE) for the three different ${\mathrm{N}\mathrm{N}\mathrm{N}-\mathrm{T}\mathrm{i}}_{S0}^{A}$ polarons (d_xz –d_yz , ${d}_{{x}^{2}-{y}^{2}}$ and d_yz ).

${\mathrm{N}\mathrm{N}\mathrm{N}-\mathrm{T}\mathrm{i}}_{S0}^{A}$	ΔE	dxy	dxz	dyz	${d}_{{x}^{2}-{y}^{2}}$	${d}_{{z}^{2}}$
dxz –dyz	0	—	24	41	—	—
${d}_{{x}^{2}-{y}^{2}}$	+25	—	10	1	50	3
dyz	+40	—	6	58	—	—

The stabilization of a different orbital symmetry can be attributed to surface effects and to the structural buckling of the TiO₂ surface structure: the tendency to stabilize the d_xz –d_yz symmetry is due to the misalignment between ${\mathrm{T}\mathrm{i}}_{S0}^{B}$ and ${\mathrm{T}\mathrm{i}}_{S0}^{A}$ atoms on the rutile surface that prevents the hybridization of the ${d}_{{x}^{2}-{y}^{2}}$ orbitals, typical for sub-surface polarons. However, surface oxygen vacancies perturb the local structure and partially suppress the buckling (see the distortions around the lattice defect in figures 1(a) and (e)): ${\mathrm{T}\mathrm{i}}_{S0}^{B}$ and ${\mathrm{T}\mathrm{i}}_{S0}^{A}$ atoms are practically aligned in the proximity of the V_O, and this alignment stabilizes a ${d}_{{x}^{2}-{y}^{2}}$ orbital symmetry for a polaron on the ${\mathrm{N}\mathrm{N}-\mathrm{T}\mathrm{i}}_{S0}^{A}$ or ${\mathrm{N}\mathrm{N}\mathrm{N}-\mathrm{T}\mathrm{i}}_{S0}^{A}$ site, in analogy with the sub-surface ${\mathrm{T}\mathrm{i}}_{S1}^{A}$ polaron.

These different orbital symmetries can have different impact on the material properties: by considering for example the surface reactivity, the DFT calculations modeling the most favorable d_xz –d_yz symmetry well reproduce the ability of polarons to hybridize with CO adsorbates and form stable complexes on the surface ${\mathrm{N}\mathrm{N}\mathrm{N}-\mathrm{T}\mathrm{i}}_{S0}^{A}$ sites [54]. These complexes are characterized by peculiar double-lobe signals in scanning tunneling microscopy (STM) measurements, well reproduced by DFT; conversely, the polaron in ${d}_{{x}^{2}-{y}^{2}}$ orbital symmetry interacts less favorably with the CO molecule, showing a weaker hybridization and less intense single-spot in the simulated STM, not supported by the experiments [65].

In the following, we focus on the dependence of the properties of the ${\mathrm{T}\mathrm{i}}_{S1}^{A}$ and the (d_xz –d_yz ) ${\mathrm{N}\mathrm{N}\mathrm{N}-\mathrm{T}\mathrm{i}}_{S0}^{A}$ polarons on different quantities, including the electronic correlation and lattice strain, and further analyze the potential impact on applications.

3.1. The U parameter

The properties of polarons depends strongly on the parameter U adopted in the DFT + U calculations to correct the electronic correlation. Figure 3 summarizes the effects of the U value on the localization of a ${\mathrm{T}\mathrm{i}}_{S1}^{A}$ polaron in terms of localized charge (evaluated by inspecting the local magnetization of the hosting site) and structural distortions. Larger U values determine a stronger localization of the polaronic charge: this is revealed by a magnetic moment progressively approaching a value of 1 μ_B (figure 3(a)). The variation in the charge localization is accompanied by a modification of the structural distortions around the polaronic sites: typically, anions (such as O atoms) around the polaron move away from localization sites as a result of the repulsive interactions with the excess electron, while cations (such as Ti atoms) are attracted by the electron polaron. Figures 3(b) and (c) show the evolution of the distance from a polaronic ${\mathrm{T}\mathrm{i}}_{S1}^{A}$ site of an octahedrally coordinated O atom (d_Ti–O) and the nearest Ti atom along [001] (d_Ti–Ti). As one might expect, the bond-length d_Ti–O increases for increasing U (figure 3(b)), due to the enhanced charge localization. Surprisingly, the attractive interaction with the nearest Ti atom is instead reduced by increasing U (d_Ti–Ti is shorter for smaller U values in figure 3(c)). This effect is shown also in figure 3(d), where we report the polaron-induced variation of the Ti–Ti distance (Δ_Ti–Ti) as compared to the Ti–Ti distance obtained in slabs with delocalized electrons: for large value of U, the polaronic distortion around the ${\mathrm{T}\mathrm{i}}_{S1}^{A}$ atom tends to the delocalized limit. This apparently counter-intuitive trend could be probably understood as a consequence of the displacement of the equatorial O oxygen atoms around the polaronic site, pushing the nearest Ti atoms away along the [001] direction.

Zoom In Zoom Out Reset image size

The (d_xz –d_yz ) ${\mathrm{N}\mathrm{N}\mathrm{N}-\mathrm{T}\mathrm{i}}_{S0}^{A}$ polarons show similar trends for the charge localization and the distortions of the coordinated O atoms (not shown), but present a different behavior as far as the nearest Ti atoms are considered. As shown in figure 3(d), the nearest Ti atoms along [001] are found at larger distance than for a delocalized electron, at variance with the attractive regime of the ${\mathrm{T}\mathrm{i}}_{S1}^{A}$ polarons. (We considered the distance between the ${\mathrm{N}\mathrm{N}\mathrm{N}-\mathrm{T}\mathrm{i}}_{S0}^{A}$ atom and the next site along [001], further away from the V_O, in order to exclude contributions from the defect.) Moreover, the extent of this polaronic distortion seems to be not particularly affected by the U parameter. These diverse characteristics are due to the d_xz –d_yz orbital symmetry of the ${\mathrm{N}\mathrm{N}\mathrm{N}-\mathrm{T}\mathrm{i}}_{S0}^{A}$ polaron: by modeling instead the surface polaron in the ${d}_{{x}^{2}-{y}^{2}}$, a similar trend for the Ti–Ti distance as obtained for the ${\mathrm{T}\mathrm{i}}_{S1}^{A}$ polaron is retrieved (not shown).

Together with the polaron localization, the energy stability is highly affected by the value of U. Interestingly, polarons on different sites are affected to a different extent. Figure 4 compares the formation energy E_POL as a function of U for ${\mathrm{N}\mathrm{N}\mathrm{N}-\mathrm{T}\mathrm{i}}_{S0}^{A}$ and ${\mathrm{T}\mathrm{i}}_{S1}^{A}$ polarons (typically, the most favorable sites for polarons formation on the rutile surface). By using the cRPA U^cRPA value, ${\mathrm{T}\mathrm{i}}_{S1}^{A}$ polarons are more favorable (by approximately 100 meV) than ${\mathrm{N}\mathrm{N}\mathrm{N}-\mathrm{T}\mathrm{i}}_{S0}^{A}$ polarons. By adopting smaller U values, the polarons become less stable and eventually delocalize (the ${\mathrm{N}\mathrm{N}\mathrm{N}-\mathrm{T}\mathrm{i}}_{S0}^{A}$ polaron is destabilized already at U ≃ 3 while the ${\mathrm{T}\mathrm{i}}_{S1}^{A}$ polaron is slightly more resistant). On the other hand, large U values favor the formation of ${\mathrm{N}\mathrm{N}\mathrm{N}-\mathrm{T}\mathrm{i}}_{S0}^{A}$ polarons: at approximately 1 eV above the U^cRPA value (U ≃ 5 eV) formation of ${\mathrm{T}\mathrm{i}}_{S1}^{A}$ polarons becomes energetically less favorable than ${\mathrm{N}\mathrm{N}\mathrm{N}-\mathrm{T}\mathrm{i}}_{S0}^{A}$. The crossing takes place at even lower U values for strongly reduced slabs (e.g., U ≃ 4.5 eV for ${c}_{{V}_{\text{O}}}\ =16.7\%$, not shown), closer to U^cRPA. Considering that the stability of the S0 polarons is not typically supported by the experimental observations (on the clean surface, at low temperature), these large U values should be avoided in the calculations. Furthermore, the curves for both polarons show a kink at U ≃ 7 eV, with E_POL following a different linear trend at higher U values: this discontinuity can be understood by inspecting the corresponding electronic properties, as discussed below.

Zoom In Zoom Out Reset image size

The E_POL trend is associated to a modification of the polaronic in-gap state. Figure 5 shows the DOS obtained for the ${\mathrm{T}\mathrm{i}}_{S1}^{A}$ polaron. The position of the polaron state can be arbitrarily controlled by varying the U: upon variations of ±1 eV around U^cRPA (up to approximately U = 5 eV), the polaron peak remains sharp and shifts within almost the whole energy gap. At higher U values, the polaronic state overlaps with the valence band, and, eventually, the d states of the hosting Ti site show no sharp features at U≳7 eV, in correspondence with the discontinuity observed for E_POL (figure 4).

Zoom In Zoom Out Reset image size

3.2. Lattice strain

The stability of small polarons depends on the strain applied to the crystal [66]. Polarons localized on ${\mathrm{T}\mathrm{i}}_{S1}^{A}$ and ${\mathrm{N}\mathrm{N}\mathrm{N}-\mathrm{T}\mathrm{i}}_{S0}^{A}$ sites are characterized by different local lattice distortions, as discussed in the previous section: while the bond-length with the coordinated O atoms are elongated by both polarons, the ${\mathrm{T}\mathrm{i}}_{S1}^{A}$ polaron tends to attract the nearest-neighbor Ti atoms, while the ${\mathrm{N}\mathrm{N}\mathrm{N}-\mathrm{T}\mathrm{i}}_{S0}^{A}$ repels the nearest Ti atoms. We can describe this scenario by considering the ${\mathrm{N}\mathrm{N}\mathrm{N}-\mathrm{T}\mathrm{i}}_{S0}^{A}$ polaron as subjected to a compressive stress, while the ${\mathrm{T}\mathrm{i}}_{S1}^{A}$ polaron experiences a tensile stress, as far as the surrounding Ti atoms are considered. This description suggests that these two inequivalent polarons could respond differently to strain. Figure 6 compares the effects of strain (applied along the [001] direction) on the energy stability E_POL for ${\mathrm{T}\mathrm{i}}_{S1}^{A}$ and ${\mathrm{N}\mathrm{N}\mathrm{N}-\mathrm{T}\mathrm{i}}_{S0}^{A}$ polarons (with respect to the low-temperature experimental lattice constant c = 2.953 Å taken as reference, i.e., 0% strain [21]). As expected, a tensile strain along [001] enhances the stability of ${\mathrm{N}\mathrm{N}\mathrm{N}-\mathrm{T}\mathrm{i}}_{S0}^{A}$ polarons to a larger extent as compared to the ${\mathrm{T}\mathrm{i}}_{S1}^{A}$ polaron. Ultimately, for very large values of the strain (≳3%), a phase transition occurs: the ${\mathrm{N}\mathrm{N}\mathrm{N}-\mathrm{T}\mathrm{i}}_{S0}^{A}$ sites become the most favorable sites for polaron formation. Even though the ${\mathrm{T}\mathrm{i}}_{S1}^{A}$ atoms remain the most favorable hosting sites up to high amount of strain ($\approx 3$%), some specific surface properties might get sensibly affected at smaller extent of strain, such as the surface reactivity (see discussion in the next section).

Zoom In Zoom Out Reset image size

We conclude this section by underling the importance of the [001] lattice parameter used in the calculations to model the TiO₂(110) surface. Due to the known tendency of GGA-based DFT + U calculations to overestimate the unit cell volume, our optimized [001] lattice parameter is about 3% larger than the experimental one, biasing the relative stability between non-equivalent polarons. By using the low-temperature experimental parameters to model the system instead of relying on the calculated lattice, the results on the polaronic properties appear more reasonably in agreement with the experimental observation.

We note also that the calculated lattice parameter depends on the reduction level of the slab: larger concentrations of oxygen vacancies show a stronger elongation of the [001] vector (3.1% for ${c}_{{V}_{\text{O}}}\ =16.7\%$) as compared to the pristine surface (3.0%). This effect is probably due to the broken bonds on the surface at the V_O sites, and might explain the tendency of rutile sample to exhibit ${\mathrm{N}\mathrm{N}\mathrm{N}-\mathrm{T}\mathrm{i}}_{S0}^{A}$ polarons at highly reducing conditions [51, 54]: the elongation along [001] mitigates the compressive stress on the ${\mathrm{N}\mathrm{N}\mathrm{N}-\mathrm{T}\mathrm{i}}_{S0}^{A}$ polaron, resulting in a stronger stability.

3.3. CO adsorption

Polarons play a decisive role for the surface reactivity of rutile TiO₂(110). The adsorption of CO molecules has been a well studied testbed case that has contributed to clarify the role of polarons in reaction processes [54, 67–74]. At low coverage, CO molecules adsorb in a vertical orientation on the TiO₂(110) surface. Adsorption on V_O sites is typically very favorable and not influenced by polarons. CO adsorption on the five fold coordinated ${\mathrm{T}\mathrm{i}}_{S0}^{A}$ is instead driven by the polarons: at low CO coverage and highly reducing conditions, the molecules preferably adsorb at ${\mathrm{N}\mathrm{N}\mathrm{N}-\mathrm{T}\mathrm{i}}_{S0}^{A}$ sites, due to the coupling with ${\mathrm{N}\mathrm{N}\mathrm{N}-\mathrm{T}\mathrm{i}}_{S0}^{A}$-polarons, as recently reported by scanning probe experiments [54]. Weakly reduced samples exhibit few CO+${\mathrm{N}\mathrm{N}\mathrm{N}-\mathrm{T}\mathrm{i}}_{S0}^{A}$-polaron complexes together with CO molecules adsorbing on sites at large distance from the V_O and from the ${\mathrm{T}\mathrm{i}}_{S1}^{A}$ polarons (${\mathrm{T}\mathrm{i}}_{S1}^{A}{\leftrightarrow}$ CO), showing no significant electronic coupling with the polaronic states.

Figure 7 shows the effects of strain on the properties of CO adsorption, as predicted by DFT calculations for a weakly reduced slab (${c}_{{V}_{\text{O}}}$ = 5.6%). As far as the low-temperature [001] lattice parameter is used in the simulations (0% of strain), DFT well describes the coexistence of CO+${\mathrm{N}\mathrm{N}\mathrm{N}-\mathrm{T}\mathrm{i}}_{S0}^{A}$-polaron complexes with ${\mathrm{T}\mathrm{i}}_{S1}^{A}{\leftrightarrow}$CO molecules (see the comparable E_ads values at 0% strain in figure 7) [54]. However, even modest tensile strain along [001] favors the polaron hopping towards the surface, and the CO+${\mathrm{N}\mathrm{N}\mathrm{N}-\mathrm{T}\mathrm{i}}_{S0}^{A}$-polaron complexes become dominant. Conversely, by applying strong compressive strain (≃−1%), the calculations predict a splitting of the CO+${\mathrm{N}\mathrm{N}\mathrm{N}-\mathrm{T}\mathrm{i}}_{S0}^{A}$-polaron complex, in favor of a polaron lying on a ${\mathrm{T}\mathrm{i}}_{S1}^{A}$ site below the CO molecule in the proximity of the V_O (labeled as CO@${\mathrm{T}\mathrm{i}}_{S1}^{A}$ in the figure) in addition to the ${\mathrm{T}\mathrm{i}}_{S1}^{A}{\leftrightarrow}$ CO configuration.

Zoom In Zoom Out Reset image size

Beside providing useful ways of controlling adsorption and catalysis processes, strain effects clearly indicate that calculations performed using the DFT + U-optimized lattice constant overestimate unrealistically the stability of CO+${\mathrm{N}\mathrm{N}\mathrm{N}-\mathrm{T}\mathrm{i}}_{S0}^{A}$-polaron complex, against experimental observation. It is clear that comparison with experiment and accurate control over the polaron properties is crucial to obtain reasonable predictions.

In summary, we reported a systematic DFT analysis on the properties of polarons on rutile TiO₂(110), focusing on the effects of lattice strain and computational details. The U parameter used in the calculation for correcting the electronic correlation on transition metals shows a broad impact, influencing the charge localization, the local lattice distortions and the electronic state of the polarons. Different types of polarons are affected to a different extent by the U value, and phase transitions could be artificially driven in the calculations depending on the adopted parameter: the sub-surface polarons typically reported by the experiments are predicted as the most stable type for U values ranging ±1 eV around the U^cRPA = 3.9 eV value calculated by cRPA simulations in the bulk. Moreover, charge trapping can be driven to different lattice sites upon applied strain: this suggests a possible mechanism for controlling the polaron localization, in addition to underline the importance of a careful parametrization in the calculations. The modeling of the excess electrons is indeed a delicate task, due to a potential energy surface characterized by numerous local minima [75], determined by the various possible polaron arrangements and the different electronic properties of these electronic states. The benchmark with experiments is of fundamental importance in order to obtain reliable predictions [76, 77] and to accurately describe complex phenomena influenced by polarons [1].

Additionally, we showed how the functionalities of polaronic materials can be tuned by controlling the polaron properties. We focused on the adsorption of CO molecules: CO adsorbates can be driven to different adsorption sites upon application of strain, arbitrarily favoring or disfavoring the coupling with surface polarons. Control over the adsorption process and the electronic properties of the adsorbates represent an interesting and coveted feature that could be exploited for applications in the field of catalysis.

Finally, we note that the development of efficient functional oxide materials requires an accurate atomic-scale understanding of the effects brought about by oxidation and reduction reactions. Our study highlights the importance of a precise modeling of reducible oxides and the need to account for different arrangement and types of polarons, which can be achieved by extensive exploration of the polaronic configuration space [27, 35, 51, 78, 79].

This work was funded by the Austrian Science Fund (FWF) project F81, SFB 'TACO'.

The data that support the findings of this study are available upon reasonable request from the authors.