Facet-dependent polaron stability in photocatalysis by SrTiO₃: a constrained DFT study

Abstract

Strontium titanate (SrTiO₃ or STO) is one of the promising photocatalysts for sustainable energy applications. Using the density functional theory (DFT) calculations, we herein study the structural and electronic factors contributing to its high photocatalytic activity and facet dependence. The constrained DFT method revealed that the hole polarons in bulk and surface STO are localized primarily on oxygen atoms. In contrast, electron polarons in bulk STO tend to delocalize over oxygen atoms unless stabilized by oxygen vacancies. The stability of hole polarons is higher at the surface O site of the (110) surface compared to the (001) surfaces. In addition, the oxygen vacancy is stable specifically at the TiO₂-terminated (001) surface. These findings provide an atomic-level insight into the relationship between polaron stability and facet dependence of photocatalysis, paving the way for the design of more efficient STO-based photocatalysts.

---

1. Introduction

Photocatalysis is one of the promising technologies for sustainable energy production, offering a potential alternative to fossil fuels and nuclear energy. For example, Al-doped SrTiO₃ (STO) photocatalyst, combined with the cocatalysts Rh/Cr₂O₃ and CoOOH, has recently been developed and demonstrated nearly 100% quantum efficiency for water splitting,^1–3 which involves breaking down water (H₂O) into hydrogen (H₂) and oxygen (O₂). This example illustrates that STO is commonly utilized in photocatalysis, not only for water splitting but also for CO₂ conversion⁴ and other applications.^5,6 While the precise design of catalysts is crucial to enhancing a photocatalytic activity, understanding the underlying reaction mechanism remains a key challenge in further improving photocatalyst performance. One of the major obstacles lies in the difficulty of accurately controlling and measuring the microscopic surface structure and the facet-dependent photocatalytic activity in experimental conditions.

In photocatalysts, the catalytic activity varies significantly across different crystal facets.^1,7,8 Controlling these facets is therefore crucial for optimizing photocatalyst performance. In the experiment of STO,¹ it was observed that a hole polaron that is employed in the oxygen-evolution reaction (OER) migrates to the (110) surface, whereas an electron polaron that is employed in hydrogen-evolution reaction (HER) migrates to the (001) surface. This highlights the significance of the stability of electron and hole polarons at the surface in understanding the photocatalytic activity of STO.

Density functional theory (DFT) calculations have been widely employed to clarify the physical properties of STO,^9,10 such as the surface structures of the catalysts.^11–16 STO has various facets, such as SrO-terminated (001) facet,¹² with varying stabilities depending on the experimental conditions.^17,18 Therefore, elucidating the facet dependence, which has been shown to be important for various photocatalysts,^7,19,20 is crucial to further improve the photocatalytic performance. However, the intricate relationship between the surface structure and the reaction mechanism—particularly the stability of polarons—remains poorly understood.

This study aims to elucidate the stability of polarons and its facet dependence to better understand the reaction mechanism behind the high photocatalytic activity in STO. Our study begins with an analysis of the bulk structure in STO, focusing on the behaviour of hole and electron polarons as well as oxygen vacancies. To simulate the localization of these polarons and the dynamics of polaron transfer, we applied the constrained DFT (CDFT) method^21–24 to charge-transfer reactions based on first-principles calculations.^7,25,26 The transmission coefficient is also calculated to account for the dynamical effects on the polaron transfer. We then explore the structure of the (001) and (110) surfaces and the facet-dependence stability of polarons to better understand the reaction mechanism of STO photocatalysts. The adsorption of a water molecule and influence of OH group at the STO surface are also considered, given the potential impact on the charge transfer.

The remainder of this paper is structured as follows: Section 2 details the computational and parameters employed for the DFT calculations. In Section 3, we present and discuss the computational results, including the bulk and surface structures, polaron stability, transmission coefficient, and adsorption of a water molecule, with a focus on understanding the facet-dependent photocatalytic activity. Finally, Section 4 summarizes the key findings of the study.

2. Computational conditions

We performed DFT calculations using the cp2k program package.²⁷ The Goedecker–Teter–Hutter pseudopotentials^28,29 and a double-ζ valence plus polarization (DZVP-MOLOPT-SR-GTH³⁰) Gaussian basis set for the orbitals were used. We employed the auxiliary augmented matrix method (ADMM),³¹ which approximates the exchange integrals to improve the speed of Hartree–Fock (HF) exchange calculations, with the ADMM basis sets (cFIT10 for Ti and Sr atoms and cFIT3 for O and H atoms). The cutoff associated with mapping gaussians onto multi-grid was 60 Ry and planewave cutoff for the finest level of the multi-grid was 400 Ry. Only the Γ point was employed for k-point sampling. The self-consistent field cycle is converged to an electronic gradient tolerance of ε_SCF = 10⁻⁵ using the orbital transformation method.³² To examine the electronic structure, the density of states (DOS) was also calculated as a function of energy which was shifted so that the Fermi level in the neutral state of STO is zero.

We employ the CDFT method in this study to localize polarons. Some previous studies have used geometric coordinates as reaction coordinates to calculate energy profiles in polaron transfer by interpolating the coordinates between reactant and product states.³³ However, this approach is limited to energetics and can hardly be applied to molecular dynamics. In contrast, the CDFT method can utilize energy-gap coordinates (defined in eqn (1) below),^25,26 which enables sampling of intermediate states in molecular dynamics simulations. The thermal fluctuations in MD simulations can influence charge transfer processes, such as polaron transfer. For example, at the TiO₂ surface, polarons can transfer among various surface sites in MD simulations at finite temperatures,³⁴ indicating the thermal fluctuations of polaron stability. The thermal fluctuations are also crucial when discussing free energies in charge-transfer reactions. For instance, in the proton transfer of aqueous silicic acid, the calculated reaction free energy is 13.9 kcal mol⁻¹—significantly lower than the enthalpy of reaction (29.3 kcal mol⁻¹) in the absence of solvent-induced fluctuations.²⁶ In addition, this capability enables the calculation of the transmission coefficient, as described below.

In the CDFT calculations, the atomic radii of the Becke density partitioning method^22,35 were 1.85 Å for Sr, 1.36 Å for Ti, 0.63 Å for O, 0.32 Å for H from ref. 36. The convergence threshold of the charge constraint is 10⁻³ electrons. The snapshots were drawn with VESTA.³⁷

2.1 Bulk

2.1.1 Cell length. The Perdew–Burke–Ernzerhof (PBE)³⁸ exchange–correlation functional was employed. To localize a polaron, which can be impossible with the generalized gradient approximation, 10.5% HF exchange that reproduces the band gap of anatase TiO₂^7,8 was included. To discuss the functional dependence on the structure, we also used the revPBE^38,39 functional and D3 dispersion correction⁴⁰ that is known to reproduce the liquid properties of water^26,41 and so on. Only in the calculations of cell length, do we employ the three (PBE, hybrid PBE, and revPBE-D3) functionals mentioned above unless otherwise noted to determine the one that best matches the experimental cell lengths. Aside from the calculations of cell length, only the hybrid PBE functional is used. It is noted that the suitability of the hybrid functional was studied in ref. 8 for TiO₂. While we limit our discussion to the three functionals mentioned above, the dependence on the choice of functional is important because there are various approaches for polaron localization, such as hybrid functionals⁸ and DFT+U method;³⁴ however, the study of functional dependence lies beyond the scope of the present work.

The initial bulk geometry of cubic perovskite STO was based on the experimental cell length of 3.8996 Å.⁴² A supercell was constructed by expanding the unit cell threefold along a, b, and c axes. The resulting supercell consists of 27 units of SrTiO₃ (135 atoms in total). Both the atomic positions and cell lengths were optimized in geometry optimizations. In the following calculations, the cell length of supercell was fixed to the optimized cell length.

2.1.2 Polaron transfer. In the simulation of a hole polaron, we added an excess positive charge (+1) to the systems, and the hole was localized on an O atom by the CDFT method. In the CDFT calculations, the constrained magnetization of an oxygen atom with the Becke density partitioning method²² was varied by starting from 0.46 to 0.52 in increments of 0.1. The value of magnetization depends on the catalysts and the sites to which a polaron is localized. For example, in TiO₂, this varied approximately from 0.4 to 0.7 at bulk sites^7,8 and from 0.7 to 0.9 at surface sites.⁷ Thus, the value of magnetization was varied from 0.46 to 0.52, because the potential energy surface exhibited a minimum in this range as discussed in Section 3.1.2. In the simulation of an electron polaron, we added an excess electron to the systems, and the excess electron was localized on a titanium atom by the CDFT method. Under these conditions, the geometry of all atoms in the system was optimized.

2.1.3 Transmission coefficient. The transmission coefficient in the hole polaron transfer was evaluated through molecular dynamics (MD) simulations of bulk STO using the computational procedure shown in Fig. 1. For this specific calculation, we employed the smaller supercell constructed by expanding the unit cell twofold along all crystal directions to reduce the computational cost of CDFT MD simulations. The normal ab initio MD (AIMD) simulation without the CDFT method was performed for 1 ps as an equilibration run in the NVT ensemble for a neutral simulation cell. In the MD simulations, the simulation time step was 2.0 fs. The time constant of thermostat^43–45 is 1 ps. A production run followed for 5 ps, during which snapshots of coordinates and velocities were stored every 1 ps, resulting in five initial conditions.

Fig. 1 Computational procedure for the equilibrium and production runs in ab initio MD (AIMD) simulations. The trajectory is denoted with a solid black line, and the initial AIMD was performed for a neutral simulation cell while the subsequent CDFT AIMD was performed for simulation cell with an excess positive charge (+1).

Subsequently, an excess positive charge (+1) was added to the systems, and CDFT MD simulations started from the above-mentioned five snapshots and performed for 1 ps in the NVT ensemble as an equilibration run. During these simulations, the magnetization constraint of 0.48 (as detailed in Section 3.1.2) was applied to a single oxygen atom. Following equilibration, 5 ps diabatic MD simulation with the CDFT method was performed, resulting in the total simulation time of 25 ps (5 initial conditions × 5 ps each).

To calculate the transmission coefficient, we employed the energy-gap coordinate. The vertical energy-gap coordinate for the CDFT method with nuclear configuration R was by^25,26

s = E_R (R) − E_P (R), (1)

in which E_R and E_P denote the total potential energies of the reactant (R) and product (P) diabatic states.

We next computed the time-correlation function (TCF) of the time derivative of the energy-gap coordinate,

C(t) = 〈ṡ(0)ṡ(t)〉, (2)

were calculated for each trajectory and averaged to obtain smooth spectra. The Fourier transform of TCF provides the dynamic analysis of vibrational couplings to the polaron transfer reaction. The method is related to the steepest-descent reaction path analysis with the projection of the coupled many degrees of freedom to the energy-gap reaction coordinate vector (see ref. 46–49 for details).

To minimize the sampling noise arising from the limited simulation time of the MD simulations, we multiplied the Gaussian window function of exp(−t²/τ²) with τ = 0.5 ps to the TCF. The TCF was truncated to 2048 steps (equivalent to 4.096 ps). To assign the vibrational modes in the spectral density, we also performed the normal mode analysis using the optimized geometry. In this calculation, the increment to be used to construct the Hessian with finite difference method was set to 0.01 bohr.

The transmission coefficient was evaluated by⁵⁰

Here, P_LZ is the Landau–Zener probability given by . Here, C is the coupling constant computed by the CDFT method and λ is the reorganization energy evaluated from the potential energy surface (see Fig. 5a). k_B is the Boltzmann constant and T is the room temperature of 298.15 K (25 °C). Ω is the effective frequency that includes the frictional effect evaluated from the spectral density (Fig. 6).^46,47,49

2.2 Surface

We examined the following surfaces that have low surface energies and high surface stability: the SrO and TiO terminations of (001) surfaces, and the O termination of (110) surface.¹² For these studies, 3D periodic boundary conditions were employed, with the cell dimensions set to 11.839 × 11.839 × 30.786 ų and 11.162 × 11.839 × 28.953 ų for the STO (001) and (110) surfaces, respectively. The slab supercells of STO surfaces were constructed by replicating a unit cell of each surface so that a polaron can be stabilized in a large supercell. The slab models of the surfaces were separated by a vacuum layer of 15 Å thickness. Along the surface normal direction, the numbers of layers are nine and eleven for the STO (001) and (110) surfaces, respectively. During geometry optimization, all atomic positions were relaxed.

Because the (001) surfaces are not stoichiometric, the non-standard formula was required to calculate the surface energy. To calculate the surface energy of the (001) surfaces, the unrelaxed cleavage energy surface was calculated as follows:¹²

where Λ represents the type of terminations, such as SrO or TiO₂. In addition, E_slab and E_bulk denote the energies of a slab and bulk STO, respectively. x, which represents the number of bulk unit cells in the SrO- and TiO₂-terminated (001) slabs, was 81. The total energy was divided by four, because the SrO- and TiO₂-terminated (001) slabs contain four surfaces in total. The unrelaxed cleavage energy disregards structural relaxation, and therefore we next calculate the relaxation energy

Finally, the surface energy was calculated by summing the unrelaxed cleavage energy and relaxation energy

E_surf(Λ) = E^(unr)_surf(Λ) + E_rel(Λ). (6)

In contrast, the surface energy at the (110) surface was computed by a standard formula

because the O-terminated (110) surface is stoichiometric. The energy difference was divided by two because the slab has two surfaces, and n_surf = 36 is the number of bulk unit cell. These surface energies were computed in the unit of eV per unit cell.

2.3 Adsorption

To calculate the adsorption strength, the adsorption energy was computed by the following definition:

E_ads = E(substrate + adsorbate) − E(substrate) − E(adsorbate), (8)

Here, E(substrate + adsorbate), E(substrate), and E(adsorbate) are the energies of the substrate and adsorbate, substrate and adsorbate, respectively. In this work, the substrate was the STO slab, and the adsorbate was a water molecule.

3. Results

3.1 Bulk

3.1.1 Cell length. Fig. 2a shows the optimized structure of neutral bulk STO, and Table 1 summarizes the optimized lattice constants (L) and band gap (E_g) using various DFT functionals. Among these, the agreement with the experimental cell length of 3.8996 Å at 296 K⁴² is the best for the hybrid PBE functional. It is noted that, if the supercell is expanded to 4 × 4 × 4 unit cell (64(SrTiO₃) with 320 atoms), the cell length changes slightly to 3.943 (hybrid PBE), 3.942 (PBE), 3.959 (revPBE-D3) Å, demonstrating the supercell's accuracy (the error is less than 0.03 Å) in this study.

Fig. 2 (a) Optimized structure of bulk STO using the hybrid PBE functional. Black lines represent the boundary of a simulation cell. (b) Projected density of states (PDOS) as a function of energy E. In panel (b), the energy axis (ordinate) is shifted such that the Fermi level in the neutral state is zero. The positive and negative DOS represent the contributions from the alpha and beta spins. Colour code: green: Sr, red: O, and light blue: Ti.

Table 1 Optimized cell length L (in Å) and band gap E_g (in eV) for various DFT functionals in STO

Functional	L	E g
a Ref. 42. b Ref. 51.
Hybrid PBE	3.946	3.12
PBE	3.968	2.79
revPBE-D3	3.960	2.81
Experiment	3.900^a	3.25 (indirect)^b
		3.75 (direct)^b

---

The experimentally measured indirect band gap energy is 3.25 eV, while the direct band gap energy is 3.75 eV.⁵¹ The agreement with the experimental band gap is also the best for the hybrid PBE functional (Table 1). These results indicate that the inclusion of HF exchange improves the agreement with the experimental physical properties. The smaller computed band gaps without HF exchange have been observed in the previous calculations using the local density approximation, where the indirect and the direct band gap energies of STO determined from the band structure calculations are 1.89 and 2.22 eV.⁵¹

Because the hybrid PBE functional more accurately reproduces the experimental physical properties, it was used for all subsequent calculations. In addition, the cell length of supercell was fixed to the optimized cell length in the following calculations.

Fig. 2b shows the projected density of states (PDOS) of bulk STO using the hybrid PBE functional. Most of the electronic states in the valence band originate from orbitals of oxygen, while most of the electronic states in the conduction band are derived from orbitals of titanium. No significant contribution from orbitals of strontium and titanium is observed near the Fermi energy.

3.1.2 Hole polaron transfer. We next set the total charge of the system to +1 for the bulk STO discussed in the previous section and applied the constraints to the magnetization of one O atom, labelled as O_D (D: donor) in Fig. 3a during geometry optimization. The value of magnetization was varied from 0.46 to 0.52. Fig. 3b shows the potential energy surface as a function of the constrained magnetization. Among the tested constraints, the constraint with magnetization of 0.48 corresponds to the lowest energy. Consequently, in the following we adopt this constraint to localize a hole polaron. The Hirshfeld spin density for the O atom is 0.50, which is close to the applied constraint value. It is worth noting that the Hirshfeld spin density and the value of constraint slightly differ because their computational methods differ. We also note that when the constraint was removed and the geometry was optimized, the Hirshfeld spin density for the O atom remained at 0.49, indicating the presence of a stable hole polaron.

Fig. 3 (a) Optimized structure of bulk STO using the hybrid PBE functional and the magnetization constraint of 0.48. A hole polaron is localized on an oxygen atom. Black lines represent the supercell of simulation. The yellow lobes represent the spin density distribution. The yellow lobes were displayed with the isosurface level of 0.0047. (b) Dependence of energy on the constrained magnetization of an O_D atom, marked in panel (a).

Fig. 3a also shows the (yellow) isosurface of spin density of the localized hole. While primarily localized on the specified oxygen atom, the hole extends also to four neighbouring oxygen atoms, each exhibiting the smaller Hirshfeld spin density of 0.10. This indicates that, in STO, a hole tends to be slightly more delocalized compared to in anatase TiO₂,⁷ in which the Hirshfeld spin density of a localized hole on an O atom typically exceeds 0.8. It is noted that if the hybrid PBE functional is changed to the PBE functional, a hole delocalizes as shown in Fig. S2a (ESI†), indicating the necessity of HF exchange for polaron localization. In addition, using the revPBE functional with 10.5% HF exchange localized a hole as shown in Fig. S2b (ESI†) and yields the magnetization of 0.51 on an O atom by Hirshfeld analysis. This shows that the different generalized gradient approximation functionals with the same fraction of HF exchange can give similar magnetizations for hole polarons, showing the robustness of a stable hole polaron in bulk SrTiO₃ within the DFT functionals employed in this work.

The distance between the donor and acceptor oxygen atoms slightly decreases from 2.79 Å in the neutral state to 2.73 Å when the polaron is localized on the donor oxygen atom to stabilize the positive charge of a polaron. On the other hand, both Sr- and Ti–O_D bond lengths respectively increase on average from 2.79 and 1.98 Å (neutral state) to 2.83 and 2.06 Å (polaron-localized state), because the positive charge of a polaron repels Sr and Ti cations. (Detailed bond lengths are provided in Fig. S1a, ESI.†) These structural adjustments promote the stabilization of the localized positive charge on the oxygen atom. Similar structural changes have also been observed in TiO₂.⁸

The PDOS in the hole-polaron-stabilized state is shown in Fig. 4. A mid-gap state associated with a hole on the O_D atom appears within the band gap, positioned slightly above the valence band maximum (VBM) at 0.69 eV. The contributions from the other atoms are close to those in Fig. 2b, indicating that the electronic structure of atoms, aside from the O_D atom, is largely unaffected by localization of a polaron.

Fig. 4 PDOS of bulk STO in the presence of a hole polaron as shown in Fig. 3a. The positive and negative DOS correspond to the alpha and beta spin states. The DOS projected on the donor O atom is magnified by ten times and marked in red filled area.

Fig. 5a shows the diabatic potential energy surface (PES) calculated via the CDFT method for the hole polaron transfer from a donor oxygen atom (O_D) to an acceptor oxygen atom (O_A). Fig. 5b and c respectively show the optimized bulk STO and the localized polaron at the transition state (TS) and product state. To calculate the diabatic potential energy surface, we adopted the following mapping potential according to the linear mixing scheme^25,26

E_α = (1 − α)E_R + αE_P. (9)

Fig. 5 (a) Diabatic (E_R and E_P) and adiabatic potential energy surface in the hole polaron transfer calculated using the CDFT method and the hybrid PBE functional. Optimized structures of the (b) TS and (c) product states. In (b), the donor (D) and acceptor (A) O atoms are also labelled. The yellow lobes represent the spin density. The yellow lobes were displayed with the isosurface levels of (b) 0.0047 and (c) 0.0050.

By varying the perturbation parameter α from 0 (reactant) to 1 (product) in constant increments of Δα, we can generate the diabatic states. The midpoint value, α = 0.5, corresponds to the TS of the reaction. The orbital in the adiabatic state was obtained by calculating the Kohn–Sham orbital using the optimized geometry through the CDFT method.

The calculated electronic coupling of polaron transfer in the bulk phase is 0.20 eV by the CDFT method. This value slightly exceeds the barrier height of 0.15 eV in the diabatic PES of Fig. 5a. In contrast, the barrier height in the adiabatic PES, also shown in Fig. 5a, is 0.05 eV. The barrier height is comparable to the thermal energy of k_BT = 0.03 eV at 298.15 K. This low activation energy is close to the experimental activation energy of 0.06–0.10 eV in N ion implanted SrTiO₃⁵² and 0.10 eV in degraded crystal of SrTiO₃.⁵³ In the calculation of the adiabatic PES, we adopted the geometry optimized with the mapping potential and obtained the Kohn–Sham orbital after removing the constraint on the magnetization. We note that quantitative comparison with experimental results is challenging for the activation energy due to differences between experiments and calculations. For instance, SrTiO₃ employed in experiments contains dopants and defects that can influence polaron transfer, while our calculations employed pristine SrTiO₃. In addition, using the linear-interpolation method,³³ the activation energy increases to 0.11 eV, which is higher than that of 0.05 eV in the adiabatic PES using the CDFT method. The lower activation energy by the CDFT method arises from the geometry relaxation at the TS.

The spin density at TS, shown in Fig. 5b, reveals that the hole is delocalized over the two oxygen atoms of the donor (D) and acceptor (A). The distance between the two O atoms of the donor and acceptor slightly increases from 2.73 Å (polaron on the donor O atom) to 2.77 Å (TS) due to repulsion between the polarons. Meanwhile, both Sr–O_D and Ti–O_D bond lengths change from 2.83 and 2.06 Å (reactant state) to 2.83 and 2.02 Å (TS) on average, respectively. These structural changes arise from the repulsion between the positively charged polaron and the Sr and Ti cations. (Individual bond lengths are provided in Fig. S1a, ESI.†)

3.1.3 Transmission coefficient. Using the calculated energy-gap coordinate in the CDFT MD simulations, we next calculated the transmission coefficient to examine the dynamical effects in polaron transfer. First, we calculated the spectral density shown in Fig. 6 for the polaron transfer in bulk STO. The spectrum mainly consists of a peak with a centre of 419 cm⁻¹. Minor contributions in the lower and higher frequencies can also be found, including small sub-peaks approximately at 208 and 688 cm⁻¹. Using the calculated spectrum, we calculated the transmission coefficient. In this calculation, the reorganization energy was estimated from the PES in Fig. 5a as λ = 0.54 eV. The calculated transmission coefficient is 1.00, mainly due to the large electronic coupling of 0.20 eV, as discussed earlier. This high value indicates the hole polaron transfer occurs via an adiabatic charge transfer mechanism.

Fig. 6 Fourier transform of the time correlation function (eqn (2)) of the time derivative of the energy-gap coordinate for the hole polaron transfer in bulk STO using the hybrid PBE functional. The ordinate is plotted as a logarithmic scale.

Normal mode analysis was conducted to assign the vibrational modes (see Fig. S3 for the illustrations, ESI†) by selecting the normal modes, the frequencies of which are the closest to the above-mentioned three peaks in Fig. 6. The lowest-frequency mode is assigned to the vibrations of the neighbouring Sr and Ti atoms, with a normal mode frequency of 205 cm⁻¹. The middle- and high-frequency modes are primarily associated with the O motions in the two neighbouring TiO₆ clusters and one TiO₆ cluster, with the normal mode frequencies of 420 and 697 cm⁻¹, respectively. These findings suggest that the collective vibrational modes of oxygen, which modulates the neighbouring OO distance and O–Ti–O angles, are strongly coupled to the polaron transfer dynamics.

3.1.4 Electron polaron. To simulate an electron polaron, we added an excess electron to the system and optimized the geometry. We initially localized an excess electron to a single Ti atom using the CDFT method with the constraint of 0.48. Afterwards, the CDFT constraint was removed, and the geometry was optimized. We found that an electron polaron could not be localized to a single Ti atom in pristine bulk STO even using the same computational procedure to localize the hole in the previous sections. This contrasts with the hole polaron that remained localized as discussed in Section 3.1.2. Furthermore, even when the HF exchange in the hybrid PBE functional is increased to 25%, which should enhance the localization of polaron,⁷ the electron polaron remains delocalized. As shown in Fig. 7a, the spin density is distributed across a single Ti layer composed of nine Ti atoms (see Fig. S4 for the views from different angles, ESI†). The average Ti–O bond length in the presence of a small polaron is 1.98 Å, identical to that in neutral pristine SrTiO₃ without an electron polaron. Additionally, the Hirshfeld charge of the Ti atoms shows minimal change, decreasing only slightly from 0.72 in the neutral state to 0.70 with an electron polaron. Such delocalization of electron polaron has also been observed both theoretically and experimentally in anatase TiO₂.^8,54 It is noted that even using the higher constraint values of 0.66 and 1.00 did not yield a localized electron in the final optimization without the CDFT constraint. Fig. 7c shows the PDOS in bulk STO with the delocalized electron polaron. In this state, a mid-gap state in the band gap slightly below the conduction band minimum (CBM) exists around at 2.6 eV. Except for the O atoms with electron polaron, the contributions from the other atoms are close to those in Fig. 2b.

Fig. 7 Optimized structure of bulk STO (a) without oxygen vacancy (V_O) and (b) with oxygen vacancy using the hybrid PBE functional. Black lines represent the boundary of a simulation cell. The yellow lobes represent the spin density and blue ones represent the spin density with the opposite spin. The yellow lobes were displayed with the isosurface levels of (a) 0.0020 and (b) 0.0035·DOS (c) without oxygen vacancy and (d) with oxygen vacancy. The DOSs in panels (c) and (d) projected on the Ti atoms on which electron polaron is delocalized are respectively magnified by three and ten times for clarity and marked in blue. TiO₅ clusters next to V_O in the (e) neutral state and the (f) electron-polaron-localized state. The bond lengths in Å are also shown, and the values in italics denote the Hirshfeld charges of Ti atoms.

We then optimized the geometry of bulk STO with one oxygen vacancy and one excess electron, because an O vacancy can induce an electron localization as observed in TiO₂.⁵⁵Fig. 7b shows the spin density of localized electrons. In contrast to the hole localized on O atoms, the electron polaron in part localized to the O vacancy site yet the other spin densities uniformly delocalize over the surrounding Ti atoms next to the oxygen vacancy. The magnetization of Ti atoms next to the oxygen vacancy ranges from −0.269 (down spin) to 0.241 (up spin). Fig. 7d shows the calculated DOS. The electron trapped state lies around 2.7 eV. Since an electron polaron does not localize in the same manner as a hole polaron and the oxygen vacancy seems to cause the localization, the stability of oxygen defects is investigated in the following surface calculations. Fig. 7e and f show the TiO₅ clusters next to O vacancy in the neutral state and the electron-polaron-localized state, respectively. Upon electron localization, the Ti–O bond lengths adjacent to O vacancy increase by 0.01–0.09 Å due to charge repulsion between the electron polaron and oxygen anions, while the Hirshfeld charge of Ti atoms next to O vacancy decreases by 0.08–0.12. This electron localization near O vacancy is facilitated by the high local structural flexibility of the lattice around O vacancy, as well as a more adaptable charge environment, which are more difficult in pristine SrTiO₃ without O vacancy as mentioned above.

3.2 Surface

3.2.1 Surface energy. Fig. 8 shows the optimized surfaces obtained from DFT calculations. At the SrO-terminated (001) surface shown in Fig. 8a, the Sr1–O1 bond length (d[Sr1–O1] = 2.80 Å) along the surface parallel direction in the surface-top layer remains close to the Sr–O bond length in bulk (2.79 Å), while d[Sr1–O2] decreases to 2.63 Å (see Fig. 8 for the atomic labels). On the other hand, d[Ti1–O1] elongates from its bulk value of 1.98 to 2.03 Å, and the O1 atom protrudes slightly toward the surface normal direction. To counterbalance the shrinkage of d[Sr1–O2], d[Sr2–O2] increases to 2.84 Å.

Fig. 8 Top and side views of the (a) SrO- and (b) TiO₂-terminated (001) and (c) STO (110) relaxed surfaces using the hybrid PBE functional. Black lines represent the simulation cell. The lower panels include the enlargement view of the surface.

The TiO₂-terminated (001) surface and the O-terminated (110) surface also exhibit the similar changes of bond lengths. At the TiO₂-terminated (001) surface in Fig. 8b, d[Sr1–O1] = 2.65 Å is also smaller than that in bulk (2.79 Å). In contrast, d[Ti1–O2] is 1.83 Å, which is smaller than that in the bulk value of 1.98 Å.

In contrast to the (001) surfaces, the (110) surface has asymmetry due to the presence of a surface-top O atom. The (110) surface has two types of surface O atoms (Fig. 8c): O1, which protrudes from the surface SrTiO layer, and O2, which is in the surface SrTiO layer. The Sr1–O1 bond length of 2.56 Å is much smaller than that in the bulk value of 2.79 Å. The (110) surface exhibits the lowest coordination number of three (two by Sr atoms and one by a Ti atom) for the surface O atoms among the three surfaces examined in this work. This is in contrast to the SrO-terminated (001) surface, where the coordination number is five (four by Sr atoms and one by a Ti atom), and the TiO₂-terminated (001) surface, where it is four (two by Sr atoms and two by Ti atoms).

The calculated surface energies of the optimized SrO- and TiO₂-(001) and (110) surfaces using the hybrid PBE functional are 1.22, 1.23, and 1.25 eV per surface unit cell. Our computational results are consistent with trends observed in the previous study using the B3PW functional,¹² where the order of the surface energy is reproduced as 1.15 eV (SrO-terminated (001) surface) <1.23 eV (TiO₂-terminated (001) surface) <1.54 eV (O-terminated (110) surface) though the differences among the facets are larger than our results, probably due to the different computational conditions. The small difference in surface energies may not be experimentally resolvable. However, both (001) and (110) surfaces have been observed in experiments,¹ indicating that these surfaces can coexist in STO. In real photocatalytic processes, which occur in aqueous environments, the stability of different facets may vary. While investigating this aspect is important, it lies beyond the scope of this study. Consequently, as all facets are plausible, this study considers all three surfaces in the subsequent discussions.

3.2.2 Adsorption energy of water. To examine the absorption of a water molecule, we calculated the adsorption energy, because the presence of hydroxyl (OH) group at surface can play a crucial role for polaron transfer, as observed in materials like TiO₂.⁷Table 2 summarizes the adsorption energies of a water molecule on the STO surfaces shown in Fig. 9. At all surfaces, the dissociative adsorption is energetically favoured, indicating that a water molecule prefers to split into H and OH species.

Table 2 Adsorption energy E_ads (in eV) for dissociative and molecular adsorptions of a water molecule computed in this study at the (001) and (110) surfaces of STO using the hybrid PBE functional

Surface	Molecular	Dissociative
SrO-(001)	−0.89	−1.44
TiO2-(001)	−0.77	−0.98
(110)	−0.70	−1.71

---

Fig. 9 Molecular and dissociative adsorptions of a water molecule adsorbed on the (a) SrO- and (b) TiO₂-terminated STO (001) and (c) (110) surfaces using the hybrid PBE functional. The enlarged views for the region enclosed by dashed lines are also shown, including the bond lengths in Å. Both molecular and dissociative adsorptions are shown. In panel (a), the definition of the depth coordinate for Fig. 10–12 is also displayed. Colour code: green: Sr, red: O, light blue: Ti, and white: H.

At the SrO-terminated (001) surface, a molecular adsorption is more stable than on the other surfaces, with an Sr–O bond length is 2.72 Å. The H atoms in the adsorbed water molecule point toward the surface normal direction. For dissociative adsorption, a water molecule splits to an H atom, which bonds to the surface O atom, and an OH group, which bridges the neighbouring surface Sr atoms with the Sr–O bond lengths of 2.55 and 2.62 Å. The computed adsorption energy for molecular adsorption closely matches the previous calculations ranging from −0.67 to −1.1 eV, while that for dissociative adsorption is smaller than the previous calculations ranging from −0.78 to −1.1 eV by various DFT functionals.¹³

For molecular adsorption at the TiO₂-terminated (001) surface, the O atom of a water molecule adsorbs to the surface Ti atom. Unlike at the SrO-terminated (001) surface, the H atoms in an adsorbed water molecule points toward the surface parallel direction. For dissociative adsorption, a water molecule splits to an OH group on the surface Ti atom and H atom on the neighbouring surface O atom. The computed adsorption energies in our DFT calculations are similar to the previous DFT calculations, ranging from −0.64 to −1.0 eV for molecular adsorption and ranging from −0.61 to −0.87 eV for dissociative adsorption using various DFT functionals.¹³

At the (110) surface, molecular adsorption can also be found on Sr atom. The adsorption energy for molecular adsorption is much larger than the dissociative case, and a water molecule strongly prefers to split on the surface. The Sr–O bond is not completely along the surface-normal direction due to the presence of the surface-top O atom. For dissociative adsorption, the dissociated proton bonds to the surface O atom, while the OH group bridges two Sr and Ti atoms. The stronger adsorption for dissociative adsorption was also reported previously.¹⁴

3.2.3 Hole polaron transfer. Fig. 10a shows the potential energy surface of a hole at the SrO and TiO₂-terminated STO (001) surfaces, as well as the STO (110) surfaces. In the calculations at surfaces, the geometry was optimized using the CDFT method and the computed energy was employed. For the SrO-terminated STO (001) surface, the energy rises to the maximum at the sub-surface O atom, and the energy difference is 0.40 eV between the bulk and surface O sites. For the TiO₂-terminated STO (001) surface, on the other hand, a hole can transfer to the sub-surface site, and also to the surface site with a stabilization energy of 0.32 eV. For the O-terminated STO (110) surface, a hole can transfer readily to the sub-surface and surface O sites with a stabilization energy of 0.37 eV. This indicates that a hole can readily transfer to the STO (110) surface. This is consistent with the experimental findings that a hole transfers to the (110) surface to induce the OER reactions.¹

Fig. 10 (a) Potential energy surface of hole polaron transfer and (b)–(d) polaron holes localized at the surface O atoms at the SrO- and TiO₂-terminated STO (001) and (110) surfaces using the hybrid PBE functional. The energy at the bulk site that is the most distant from the surface was set to zero. The definition of the depth coordinate is displayed in Fig. 9a. The yellow lobes were displayed with the isosurface levels of (b) 0.0030, (c) 0.0033, and (d) 0.0031.

As shown in Fig. 10b, when a hole is at the outermost O atom of the SrO-terminated STO (001) surface, the Ti–O bond length slightly decreases by 8% as the Ti atom moves toward the surface, while Sr–O bond length remains almost the same. For the TiO₂-terminated STO (001) surface as shown in Fig. 10c, Ti–O bond length increases by 5%, and the protrusion of O atom along the surface normal direction is pronounced. When a hole is localized on the surface O atom at the STO (110) surface, Ti- and Sr–O bond lengths respectively increase from 1.74 and 2.56 Å to 1.84 (6% increase) and 2.73 Å (7% increase), as shown in Fig. 10d. Such structural changes in terms of bond lengths upon localization of a hole were also observed in bulk discussed in Section 3.1.2. The most stable polaron at the STO (110) surface may be related with the flexibility of surface O atom. In fact, the surface O atom of the STO (110) surface has the lowest coordination number of three as discussed in Section 3.2.1, and therefore is thought to be the most flexible among the three surfaces considered in this study.

To investigate the impact of hydroxylation, we computed the potential energy profile as a function of polaron position at the STO surfaces on which one water molecule split into OH and H (see the Section 3.2.2 for the adsorption structures). The energy profile in Fig. 11 remains largely unchanged from that at the pristine surfaces, except at the surface OH site. At the TiO₂-terminated (001) surface, however, the energy at the surface oxygen site decreases compared to the pristine surface. However, the polaron transfer from the surface O site to the OH group is energetically unfavourable with an increase of 0.13 eV. At the SrO-terminated (001) surface, a polaron on the OH group is the most stable among the OH groups at the three surfaces studied in this work. Nevertheless, the polaron transfer should be highly unlikely due to the instability of polaron at the surface O site. In the case of the (110) surface, the energy rises to 0.14 eV at the OH site. This suggests that the hole transfer proceeds through the surface oxygen site rather than an OH group during oxidation reactions. This contrasts with the case of anatase TiO₂,⁷ where the presence of OH group stabilizes a hole polaron.

Fig. 11 Same as Fig. 10a, yet OH group is attached the surfaces from water adsorption using the hybrid PBE functional.

We further computed the potential energy surface for various positions of oxygen vacancies at the STO surfaces, as shown in Fig. 12. At all surfaces, the energy drops at the surface O sites. This would be because the structural changes at surface are easier than in bulk. The TiO₂-terminated (001) surface exhibits the more pronounced energy drop than the other surfaces. This indicates that O vacancy can be formed more readily at the TiO₂-terminated (001) surface. A higher concentration of O vacancies at the TiO₂-terminated (001) surface implies more trapping sites for electron polaron, consistent with the experimental findings that an electron polaron transfers to the (001) surface to induce the HER reactions.¹

Fig. 12 Potential energy surface for various positions of oxygen vacancies at the STO surfaces using the hybrid PBE functional.

4. Conclusions

This study provides insights into the photocatalytic activity of SrTiO₃ (STO), particularly focusing on its facet-dependent performance in photocatalysis. By employing density functional theory (DFT) calculations, including constrained DFT (CDFT) methods, we simulated the localization and transfer of hole and electron polarons in bulk STO and on its various surfaces.

Our investigation into the bulk STO structure revealed that the hybrid PBE functional closely more accurately reproduces experimental physical properties of STO, particularly in terms of cell length and band gap energy, making it a reliable choice for subsequent calculations. The simulation of hole polaron transfer demonstrated that the process in bulk STO has a low activation barrier, with a calculated barrier height close to the thermal energy at room temperature. In bulk STO, a hole polaron localizes on O atoms, while an electron polaron can localize with the existence of an oxygen vacancy. The dynamical effects of polaron transfer were also examined through by computing the transmission coefficient, and the hole polaron transfer occurs through an adiabatic charge transfer mechanism.

When exploring the surface properties, we considered both SrO and TiO₂-terminated (001) surfaces and the O-terminated (110) surface, which are known for their stability. The calculations of surface energies and adsorption energies provided a quantitative understanding of the surface stability and the interaction with a water molecule and the STO surface. The dissociative adsorption of a water molecule is energetically favoured, and the computed adsorption energies are consistent with the previous DFT studies. At the STO surfaces, a hole is localized on the surface O atom of the (110) surface most stably, and an OH group the surface does not enhance the stability. In contrast, the oxygen vacancy is most stable on the TiO₂-terminated (001) surface. The insights gained from this study can inform future efforts in developing advanced STO photocatalysts for sustainable energy applications.

Author contributions

This manuscript was written through the contributions of all authors. All authors have given approval to the final version of the manuscript.

Data availability

The authors confirm that the data supporting the findings of this study are available within the article and/or its ESI.†

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This study was financially supported by the Iketani Science and Technology Foundation. We used supercomputers from ACCMS, Kyoto University.

References

1. T. Takata, J. Jiang, Y. Sakata, M. Nakabayashi, N. Shibata, V. Nandal, K. Seki, T. Hisatomi and K. Domen, Nature, 2020, 581, 411–414.
2. H. Nishiyama, T. Yamada, M. Nakabayashi, Y. Maehara, M. Yamaguchi, Y. Kuromiya, Y. Nagatsuma, H. Tokudome, S. Akiyama, T. Watanabe, R. Narushima, S. Okunaka, N. Shibata, T. Takata, T. Hisatomi and K. Domen, Nature, 2021, 598, 304–307.
3. Y. Ham, T. Hisatomi, Y. Goto, Y. Moriya, Y. Sakata, A. Yamakata, J. Kubota and K. Domen, J. Mater. Chem. A, 2016, 4, 3027–3033.
4. T. Nakamoto, S. Iguchi, S. Naniwa, T. Tanaka and K. Teramura, Catal. Sci. Technol., 2023, 13, 4534–4541.
5. K. Nakashima, H. Takahama, M. Yoshida, K. Yamaguchi and K. Hata, Inorg. Chem., 2024, 63, 44–49.
6. G. Zvejnieks, L. L. Rusevich, D. Gryaznov and E. A. Kotomin, Phys. Chem. Chem. Phys., 2019, 21, 23541–23551.
7. T. Joutsuka, H. Yoshinari and S. Yamauchi, Bull. Chem. Soc. Jpn., 2021, 94, 106–111.
8. A. R. Elmaslmane, M. B. Watkins and K. P. McKenna, J. Chem. Theory Comput., 2018, 14, 3740–3751.
9. A. Janotti, J. B. Varley, M. Choi and C. G. Van de Walle, Phys. Rev. B:Condens. Matter Mater. Phys., 2014, 90, 085202.
10. D. A. Tenne, I. E. Gonenli, A. Soukiassian, D. G. Schlom, S. M. Nakhmanson, K. M. Rabe and X. X. Xi, Phys. Rev. B:Condens. Matter Mater. Phys., 2007, 76, 024303.
11. E. Heifets, S. Piskunov, E. A. Kotomin, Y. F. Zhukovskii and D. E. Ellis, Phys. Rev. B:Condens. Matter Mater. Phys., 2007, 75, 115417.
12. R. I. Eglitis and D. Vanderbilt, Phys. Rev. B:Condens. Matter Mater. Phys., 2008, 77, 195408.
13. E. Holmström and A. S. Foster, J. Chem. Theory Comput., 2017, 13, 6301–6307.
14. Z. Wang, X. Hao, S. Gerhold, Z. Novotny, C. Franchini, E. McDermott, K. Schulte, M. Schmid and U. Diebold, J. Phys. Chem. C, 2013, 117, 26060–26069.
15. R. Garcia-Diaz, M. T. Romero de la Cruz, R. Ochoa Valiente, J. Guerrero-Sanchez and G. Hernández Cocoletzi, Appl. Surf. Sci., 2019, 487, 1394–1402.
16. M. Sokolov, Y. A. Mastrikov, G. Zvejnieks, D. Bocharov, V. Krasnenko, K. S. Exner and E. A. Kotomin, Adv. Theory Simul., 2023, 6, 2200619.
17. J. Ciston, H. G. Brown, A. J. D’Alfonso, P. Koirala, C. Ophus, Y. Lin, Y. Suzuki, H. Inada, Y. Zhu, L. J. Allen and L. D. Marks, Nat. Commun., 2015, 6, 7358.
18. S. Kawasaki, E. Holmström, R. Takahashi, P. Spijker, A. S. Foster, H. Onishi and M. Lippmaa, J. Phys. Chem. C, 2017, 121, 2268–2275.
19. K.-I. Katsumata, Y. Okayasu, N. Matsushita and K. Okada, Chem. Lett., 2013, 42, 618–620.
20. C. Y. Toe, M. Lamers, T. Dittrich, H. A. Tahini, S. C. Smith, J. Scott, R. Amal, R. van de Krol, F. F. Abdi and Y. H. Ng, Mater. Adv., 2022, 3, 2200–2212.
21. Q. Wu and T. Van Voorhis, Phys. Rev. A:At., Mol., Opt. Phys., 2005, 72, 024502.
22. N. Holmberg and K. Laasonen, J. Chem. Theory Comput., 2017, 13, 587–601.
23. P. H. Dederichs, S. Blügel, R. Zeller and H. Akai, Phys. Rev. Lett., 1984, 53, 2512–2515.
24. B. Kaduk, T. Kowalczyk and T. Van Voorhis, Chem. Rev., 2012, 112, 321–370.
25. T. Joutsuka and K. Ando, Chem. Lett., 2021, 50, 1325–1328.
26. T. Joutsuka and K. Ando, J. Phys. Chem. B, 2020, 124, 8323–8330.
27. T. D. Kühne, M. Iannuzzi, M. D. Ben, V. V. Rybkin, P. Seewald, F. Stein, T. Laino, R. Z. Khaliullin, O. Schütt, F. Schiffmann, D. Golze, J. Wilhelm, S. Chulkov, M. H. Bani-Hashemian, V. Weber, U. Borštnik, M. Taillefumier, A. S. Jakobovits, A. Lazzaro, H. Pabst, T. Müller, R. Schade, M. Guidon, S. Andermatt, N. Holmberg, G. K. Schenter, A. Hehn, A. Bussy, F. Belleflamme, G. Tabacchi, A. Glöβ, M. Lass, I. Bethune, C. J. Mundy, C. Plessl, M. Watkins, J. VandeVondele, M. Krack and J. Hutter, J. Chem. Phys., 2020, 152, 194103.
28. S. Goedecker, M. Teter and J. Hutter, Phys. Rev. B:Condens. Matter Mater. Phys., 1996, 54, 1703–1710.
29. C. Hartwigsen, S. Goedecker and J. Hutter, Phys. Rev. B:Condens. Matter Mater. Phys., 1998, 58, 3641–3662.
30. J. VandeVondele and J. Hutter, J. Chem. Phys., 2007, 127, 114105.
31. M. Guidon, J. Hutter and J. VandeVondele, J. Chem. Theory Comput., 2010, 6, 2348–2364.
32. J. VandeVondele and J. Hutter, J. Chem. Phys., 2003, 118, 4365–4369.
33. N. A. Deskins and M. Dupuis, J. Phys. Chem. C, 2009, 113, 346–358.
34. P. M. Kowalski, M. F. Camellone, N. N. Nair, B. Meyer and D. Marx, Phys. Rev. Lett., 2010, 105, 146405.
35. A. D. Becke, J. Chem. Phys., 1988, 88, 2547–2553.
36. P. Pyykkö and M. Atsumi, Chem. – Eur. J., 2009, 15, 186–197.
37. K. Momma and F. Izumi, J. Appl. Crystallogr., 2011, 44, 1272–1276.
38. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
39. Y. Zhang and W. Yang, Phys. Rev. Lett., 1998, 80, 890.
40. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
41. T. Joutsuka, J. Phys. Chem. B, 2022, 126, 4565–4571.
42. Y. A. Abramov, V. G. Tsirelson, V. E. Zavodnik, S. A. Ivanov and I. D. Brown, Acta Crystallogr., Sect. B:Struct. Sci., 1995, 51, 942–951.
43. S. Nosé, J. Chem. Phys., 1984, 81, 511–519.
44. W. G. Hoover, Phys. Rev. A:At., Mol., Opt. Phys., 1985, 31, 1695–1697.
45. G. J. Martyna, M. L. Klein and M. Tuckerman, J. Chem. Phys., 1992, 97, 2635–2643.
46. K. Ando and S. Kato, J. Chem. Phys., 1991, 95, 5966–5982.
47. K. Ando, J. Chem. Phys., 1994, 101, 2850–2862.
48. T. Joutsuka and K. Ando, J. Phys. Chem. A, 2011, 115, 671–677.
49. T. Joutsuka and K. Ando, J. Phys. Chem. A, 2011, 115, 678–684.
50. J. E. Straub and B. J. Berne, J. Chem. Phys., 1987, 87, 6111–6116.
51. K. van Benthem, C. Elsässer and R. H. French, J. Appl. Phys., 2001, 90, 6156–6164.
52. A. Bhogra, A. Masarrat, R. Meena, D. Hasina, M. Bala, C.-L. Dong, C.-L. Chen, T. Som, A. Kumar and A. Kandasami, Sci. Rep., 2019, 9, 14486.
53. W. Liu, G.-Y. Yang and C. A. Randall, Jpn. J. Appl. Phys., 2009, 48, 051404.
54. M. Setvin, C. Franchini, X. Hao, M. Schmid, A. Janotti, M. Kaltak, C. G. Van de Walle, G. Kresse and U. Diebold, Phys. Rev. Lett., 2014, 113, 086402.
55. A. J. Tanner and G. Thornton, J. Phys. Chem. Lett., 2022, 13, 559–566.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d4cp04725g

---