Ab initio description of oxygen vacancies in epitaxially strained \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {SrTiO}_{{3}}$$\end{document}SrTiO3 at finite temperatures

Epitaxially grown \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {SrTiO}_{{3}}$$\end{document}SrTiO3 (STO) thin films are material enablers for a number of critical energy-conversion and information-storage technologies like electrochemical electrode coatings, solid oxide fuel cells and random access memories. Oxygen vacancies (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathrm{V}_{{\mathrm{O}}}}$$\end{document}VO), on the other hand, are key defects to understand and tailor many of the unique functionalities realized in oxide perovskite thin films. Here, we present a comprehensive and technically sound ab initio description of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathrm{V}_{{\mathrm{O}}}}$$\end{document}VO in epitaxially strained (001) STO thin films. The novelty of our first-principles study lies in the incorporation of lattice thermal excitations on the formation energy and diffusion properties of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathrm{V}_{{\mathrm{O}}}}$$\end{document}VO over wide epitaxial strain conditions (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-4 \le \eta \le +4$$\end{document}-4≤η≤+4%). We found that thermal lattice excitations are necessary to obtain a satisfactory agreement between first-principles calculations and the available experimental data for the formation energy of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathrm{V}_{{\mathrm{O}}}}$$\end{document}VO. Furthermore, it is shown that thermal lattice excitations noticeably affect the energy barriers for oxygen ion diffusion, which strongly depend on \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta $$\end{document}η and are significantly reduced (increased) under tensile (compressive) strain. The present work demonstrates that for a realistic theoretical description of oxygen vacancies in STO thin films is necessary to consider lattice thermal excitations, thus going beyond standard zero-temperature ab initio approaches.

(DFT) due to Perdew, Burke, and Ernzerhof (GGA-PBE) 42 as is implemented in the VASP software 43 . A "Hubbard-U" scheme 44 was employed for a better treatment of the Ti 3d electronic orbitals with a selected U value of 2.0 eV. (The main conclusions presented in this article do not appreciably depend on this particular choice as demonstrated by numerical tests carried out for U = 4 eV, see Supplementary Fig. 1 and "Formation energy of oxygen vacancies at T = 0" section.) The energy band gap of (001) SrTiO 3 thin films was accurately estimated with the range-separated hybrid functional HSE06 45 for the equilibrium structures previously determined at the GGA-PBE+U level. We used the "projector augmented wave" method 46 to represent the ionic cores and considered the following electronic states as valence: Sr 4s, 4p and 5s; Ti 3p, 4s and 3d; O 2s and 2p. Wave functions were represented in a plane-wave basis truncated at 650 eV.
For simulation of the stoichiometric systems, we employed a 20-atoms simulation cell that allows to reproduce the usual ferroelectric and O 6 antiferrodistortive (AFD) distortions in perovskite oxides 8,16 (Fig. 1). Off-stoichiometric systems containing oxygen vacancies, V O , were generated by removing oxygen atoms from either equatorial (Eq) or apical (Ap) positions (Fig. 1). Simulation cells of different sizes were considered in order to quantify the effects of oxygen vacancy concentration on the obtained formation energy results. In particular, the following compositions were investigated: Sr 4 ( SrTiO 2.94 ). For integrations within the Brillouin zone (BZ), we used a Ŵ-centered k-point grid of 6 × 8 × 8 for the 20-atoms simulation cell and scaled it conveniently to maintain an equivalent k-point density for the rest of cases. All oxygen vacancies were assumed to be neutrally charged and non-magnetic ( V O ) as in previous DFT studies this configuration has been shown to render the lowest energy for bulk off-stoichiometric SrTiO 3 9,47,48 .
The geometry relaxations of epitaxially strained (001) SrTiO 3 and SrTiO 3−δ were carried out by using a conjugated gradient algorithm that allowed to change the simulation-cell volume and atomic positions while constraining the length and orientation of the two in-plane lattice vectors (that is, |a| = |b| and γ = 90 • ). Periodic boundary conditions were applied along the three lattice-vector directions, thus the influence of surface effects were systematically neglected in our simulations. This type of calculations are known as "strained-bulk" geometry relaxations and typically are considered to be a good approximation for thin films presenting thicknesses of at least few nanometers 8,[15][16][17]19 . The simulated systems were assumed to be elastically coupled to a substrate thus the existence of possible stress relaxation mechanisms in the thin films were also neglected. The geometry relaxations were stopped when the forces on the ions were smaller than 0.01 eV/Å. By using these parameters we obtained zero-temperature energies that were converged to within 0.5 meV per formula unit.
The electric polarization of stoichiometric and off-stoichiometric (001) SrTiO 3 thin films were estimated with the Born effective charges method 8,16 . In this approach, the electric polarization is calculated via the formula: where is the volume of the cell, κ runs over all the atoms, α, β = x, y, z represent the Cartesian directions, u κ is the displacement vector of the κ-th atom as referred to a non-polar reference phase, and Z * κ the Born effective charge tensor calculated for a non-polar reference state. It is worth noting that the presence of oxygen vacancies typically induced a notable reduction in the energy band gap of off-stoichiometric systems, which in some cases led to the appearence of metallic states. Consequently, estimation of the electric polarization with the more accomplished and accurate Berry phase formalism was not possible for all the analyzed compositions and thus we opted for systematically using the approximate Born effective charges method 8,49 . Phonon calculations. To estimate phonon frequencies we employed the "small-displacement" approach 50 , in which the force-constant matrix of the crystal is calculated in real space by considering the proportionality between the atomic displacements and forces when the former are sufficiently small (in the present study this condition was satisfied for atomic displacements of 0.02 Å). Large supercells containing 160 atoms were employed to guarantee that the elements of the force-constant matrix presented practically negligible values at the largest atomic separations. We used a dense k-point grid of 3 × 3 × 3 for the calculation of the atomic forces with VASP. The computation of the nonlocal parts of the pseudopotential contributions were performed in reciprocal space in order to maximise the numerical accuracy. Once a force-constant matrix was determined, we Fourier transformed it to obtain the phonon frequencies for any arbitrary k-point in the first BZ. This latter step was performed with the PHON code 50 , in which the translational invariance of the system is exploited to ensure that the three acoustic branches are exactly zero at the Ŵ point. Central differences for the atomic forces, that is, both positive and negative atomic displacements, were considered. A complete phonon calculation involved the evaluation of atomic forces for 120 (114) different stoichiometric (off-stoichiometric) configurations with the technical parameters just described. In order to accurately compute F qh vac (see below), we employed a dense k -point grid of 16 × 16 × 16 for BZ integration. With these settings we found that the calculated quasi-harmonic free energies were accurate to within 5 meV per formula unit. Examples of full phonon spectra calculated for epitaxially strained stoichiometric and non-stoichiometric STO are provided in the Supplementary Fig. 2. Free energy calculations. We computed the quasi-harmonic Gibbs free energy associated with the formation of neutral oxygen vacancies, G qh vac , as a function of epitaxial strain, η ≡ (a − a 0 )/a 0 (where a 0 represents the equilibrium in-plane lattice parameter calculated for the stoichiometric system), and temperature, T, by following the approach introduced in previous works 25,41 . Next, we briefly summarize the key aspects and technical details of the employed quasi-harmonic Gibbs free energy calculation method.
(1)  www.nature.com/scientificreports/ The formation Gibbs free energy of non-magnetic and neutrally charged V O can be expressed as 25,41 : where subscript "vac" indicates the quantity difference between the off-stoichiometric and stoichiometric systems (e.g., E vac ≡ E SrTiO 3−δ − E SrTiO 3 ), E vac accounts for the static contributions to the free energy (i.e., calculated at T = 0 conditions by considering the atoms fixed at their equilibrium lattice positions 8 ), F qh vac for the vibrational contributions to the free energy, and µ O is the chemical potential of free oxygen atoms. The vibrational free energy of stoichiometric and off-stoichiometric systems were estimated with the quasi-harmonic formula 51-55 : where N q is the total number of wave vectors used for integration within the Brillouin zone and the dependence of the phonon frequencies, ω qs , on epitaxial strain is explicitly noted.
It is well known that first-principles estimation of µ O with DFT+U methods is challenging and may lead to large errors 56,57 . Such inherent limitations make the prediction of V O formation energies by exclusively using DFT approaches difficult and probably also imprecise. Notwithstanding, since (i) the chemical potential of free oxygen atoms does not depend on the epitaxial strain conditions but on the experimental thin film synthesis conditions, and (ii) our main goal here is to unravel the impact of lattice excitations on the formation energy and diffusion of V O as a function of η , we can safely base our analysis on the results obtained for the thermodynamically shifted Gibbs free energy: In other words, rather than adopting experimental values for µ O and/or applying empirical corrections to the calculated vacancy formation energies 40,57 , here we select arbitrary values for the oxygen-gas chemical potential without any loss of generality.
Nudged elastic band calculations. Ab initio nudged-elastic band (NEB) calculations 58 were performed to estimate the activation energy for V O diffusion in expitaxially strained (001) SrTiO 3 at zero temperature. Our NEB calculations were performed for reasonably large 2 × 2 × 2 or 3 × 3 × 3 supercells containing several tens of atoms 47 . We used q-point grids of 8 × 8 × 8 or 6 × 6 × 6 and an energy plane-wave cut-off of 650 eV. Six intermediate images were used to determine the energy barrier for oxygen diffusion along the most likely V O diffusion paths in the absence of thermal excitations. The geometry optimizations were halted when the total forces on the atoms were smaller than 0.01 eV · Å −1 . The NEB calculations were performed for five epitaxialstrain equidistant points in the interval −4 ≤ η ≤ 4%.
Ab initio molecular dynamics simulations. First-principles molecular dynamics (AIMD) simulations based on DFT were performed in the canonical (N, V, T) ensemble. The selected volumes and geometries were those determined at zero-temperature conditions, hence we neglected thermal expansion effects. The concentration of oxygen vacancies in the off-stoichiometric systems was also considered to be independent of T and equal to ≈ 1.6 %. The temperature in the AIMD simulations was kept fluctuating around a set-point value by using Nose-Hoover thermostats. Large simulation boxes containing 317 atoms ( Sr 64 Ti 64 O 189 ) were employed in all the AIMD simulations and periodic boundary conditions were applied along the three Cartesian directions. Newton's equation of motion were integrated by using the customary Verlet's algorithm and a time-step length of δt = 10 −3 ps. Ŵ-point sampling for integration within the first Brillouin zone was employed in all the AIMD simulations. The calculations comprised total simulation times of t total ∼ 10 ps. We performed three AIMD simulations at T = 1000 , 1500, and 2000 K for off-stoichiometric STO thin films considering epitaxial strains of −3.6 , 0 and +3.6%.
The mean square displacement (MSD) of oxygen ions was estimated with the formula 59 : where r i (t j ) is the position of a migrating ion i at time t j ( = j · δt ), τ represents a lag time, n τ = τ/δt , N ion is the total number of mobile ions, and N step the total number of time steps. The maximum n τ was chosen equal to N step /3 (i.e., equivalent to ∼ 3-4 ps), hence we could accumulate enough statistics to reduce significantly the MSD(τ ) fluctuations at the largest τ (see the error bars in the MSD plots presented in the following sections). Oxygen diffusion coefficients were subsequently obtained with the Einstein relation: The T-dependence of the oxygen diffusion coefficient was assumed to follow the Arrhenius formula: 6τ . www.nature.com/scientificreports/ where D 0 is known as the pre-exponential factor, E a is the activation energy for ionic diffusion, and k B the Boltzmann constant.

Results and discussion
We start by discussing the zero-temperature phase diagram of stoichiometric (001) STO thin films calculated with first-principles methods. The changes in the structural and electric polarization properties induced by the presence of equatorial (Eq) and apical (Ap) oxygen vacancies ( V O ) are subsequently explained. The impact of thermal effects on the formation energy of V O is analyzed for a wide range of epitaxial strain ( −4 η +4 %) and temperature ( 0 ≤ T ≤ 1000 K) conditions. We also report and compare the energy barriers for ionic oxygen diffusion in (001) STO thin films estimated by neglecting and by taking into account T-induced lattice vibrations. Insightful connections between our theoretical results and experimental measurements are provided whenever the latter are available.
It is worth noting that in experiments the epitaxial strains that are usually achieved are of the order of 1-2% since under larger |η| 's oxide perovskites tend to elastically relax via the formation of dislocations and other structural imperfections 60 . Nevertheless, misfit strain relaxation in oxide perovskites can be precluded, to some extent, through modulation of the thin film thickness and compositional engineering 61,62 . For instance, an experimentally attained tensile epitaxial strain of +3.9 % has been recently reported for oxygen-deficient SrCoO 3−δ thin films 25 . The epitaxial strain conditions considered in the present computational study, therefore, aim to reproduce the broad |η| ranges realized in experiments. Figure 2 shows the structural, electric polarization and energy band gap properties of stoichiometric (001) SrTiO 3 thin films estimated with first-principles methods (i.e., density functional theory -DFT-, "Computational methods" section) at zero temperature. Different crystalline phases are stabilized as a result of varying the epitaxial strain conditions, which are described in detail next. We note that several authors have previously reported analogous DFT results to ours 15,[63][64][65] and that the best agreement with the present calculations is obtained for work 65 , in which antiferrodistortive oxygen octahedra rotations (AFD) were also explicitly modeled.

Zero-temperature properties of stoichiometric (001) SrTiO 3 thin films.
In the epitaxial strain interval η −2 %, we observe the stabilization of a tetragonal I4cm phase that is characterized by a significant out-of-plane polarization ( P z , Fig. 2a) and antiphase out-of-plane O 6 rotations ( AFD − z , Fig. 2c). Coexistence of the order parameters P z and AFD − z is quite unique as they normally tend to oppose each other 31 , a polar-antiferrodistortive interplay that has been experimentally observed and characterized as a function temperature and η 30 . Under large compressive strain, half of the Ti-O bond lengths involving oxygen atoms in apical positions are significantly elongated as compared to those involving O ions in equatorial positions (Fig. 2b), a structural distortion that signals the presence of out-of-plane polarization 16,18 . In the epitaxial strain interval −2 η 0 %, a tetragonal I4/mcm phase appears that presents null electric polarization (Fig. 2a) and moderate antiphase out-of-plane O 6 rotations (Fig. 2c). In this phase, the length of the Ti-O bonds are all pretty similar regardless of the positions occupied by the oxygen atoms (Fig. 2b). It is worth noting that when some tiny monoclinic lattice distortions in the generated equilibrium geometries (i.e., α ∼ 0.1 degrees) are not disregarded the identification of this phase is also compatible with a non-polar C2/c phase that has been recently predicted for metallic LaNiO 3 thin films 66 .
In the epitaxial strain interval 0 η +2 %, a noticeable in-plane electric polarization, P xy , appears in the system that coexists with small AFD − z O 6 rotations (Fig. 2a,c). The resulting crystal phase is orthorhombic and its symmetry can be ascribed to the polar space group Ima2. Under tensile strain, half of the Ti-O bond lengths involving oxygen atoms in equatorial positions are significantly elongated as compared to those involving O ions in apical positions (Fig. 2b), a structural distortion that produces a significant in-plane polarization 16,18 . In the epitaxial strain interval η +2 %, the antiphase out-of-plane O 6 rotations completely disappear and P xy grows steadily under increasing epitaxial strain. In this latter case, the optimized crystal structure is also orthorhombic and its symmetry can be identified with the space group Amm2. Figure 2d shows the energy band gap of (001) SrTiO 3 thin films, E gap , estimated as a function of epitaxial strain with the range-separated hybrid functional HSE06 45 . The reason for including this information here will become clearer in the next subsection, where we explain the oxygen vacancy formation energy results obtained at zero temperature. It is worth noting that E gap increases noticeably under either tensile or compressive biaxial strain as compared to the corresponding zero-strain value. For instance, at η = 0 the energy band gap amounts to 3.2 eV whereas at η = ±4 % is approximately equal to 3.9 eV. Such a η-induced E gap trend is markedly different from the one predicted for binary oxides like CeO 2 and TiO 2 by using analogous first-principles methods 21 , which displays a significant E gap reduction under tensile biaxial strain. The reason for such a difference in E gap behaviour is likely to be related to the larger changes in the dielectric susceptibility that can be induced by epitaxial strain in STO thin films as compared to binary oxides 21,63 . Formation energy of oxygen vacancies at T = 0. Figure 3 shows the formation energy of oxygen vacancies calculated for (001) STO thin films at zero temperature, E vac . The concentration of V O considered in this case renders the composition SrTiO 2.75 (analogous vacancy energy results obtained for smaller oxygen vacancy concentrations are explained below). A small decrease in E vac is observed as the biaxial strain changes from compressive to tensile (i.e., of > 1 % when considering the two limiting cases η = ±4 %, Fig. 3a). For most www.nature.com/scientificreports/ η cases, it seems more favourable to create apical V O than equatorial, however the formation energy differences between the two cases are pretty small (i.e., | E vac | ≈ 0.01 eV per formula unit, Fig. 3a). We repeated the E vac (η) calculations shown in Fig. 3a by adopting a larger value of the technical parameter U (i.e., equal to 4 eV, "General technical details" section) and found quantitatively equivalent results to the ones just explained for U = 2 eV (Supplementary Fig. 1). In particular, the creation of apical oxygen vacancies was found to be energetically most favourable for most of the investigated η 's and the formation energy differences between apical and equatorial V O 's were also of the order of 0.01 eV. Moreover, analogous E vac results were obtained by (i) considering the presence of long-range dispersion interactions [67][68][69][70] and (ii) employing other exchange-correlation DFT functionals like the popular local density approximation (LDA 71 ) and the recently proposed meta-GGA functional SCAN 72,73 (Supplementary Fig. 1). The study of the η-dependence of the formation energy of oxygen vacancies in STO thin films, therefore, appears to be quite insensitive to the particular choice of the U parameter and even also of the adopted DFT exchange-correlation functional. This outcome is likely to be a consequence of large error cancellations between the energies calculated for stoichiometric and non-stoichiometric systems 8 . Hereafter, all the reported results correspond to PBE+U ( U = 2 eV) calculations.
For negative η values, the creation of oxygen vacancies in either equatorial (Fig. 3c) or apical (Fig. 3d) positions has a dramatic effect on the electric polarization of the system. In particular, the sizable out-of-plane polarization found in stoichiometric STO thin films (Fig. 2a) practically disappears when V O are exclusively created in apical positions. Meanwhile, when oxygen vacancies are generated solely in equatorial positions a non-negligible inplane polarization of ≈ 7 µC cm −2 appears for any value of compressive epitaxial strain. For positive η values, on the other hand, the general behaviour of the electrical polarization is quite similar to that found for the analogous stoichiometric thin films, although the size of P xy appreciably decreases ( ∼ 10%). Figure 3b shows the volume difference between SrTiO 2.75 and stoichiometric thin films, V vac , expressed as a function of epitaxial strain. The creation of neutral oxygen vacancies in oxide perovskites typically induces an increase in volume, the so-called chemical expansion, due to the electronic reduction of transition metal ions that are located close to V O 's 41,74 . For present purposes, it is interesting to analyze the η-dependence of V vac www.nature.com/scientificreports/ because this quantity has been found to be correlated with the contribution of lattice thermal excitations to the formation energy of V O at finite temperatures 41 . As regards equatorial oxygen vacancies, V vac turns out to be positive and moderately large (small) under tensile (compressive) epitaxial strain. By contrast, the creation of apical oxygen vacancies is accompanied by negative (positive) and large V vac absolute values (small) at large compressive (tensile) epitaxial strain (Fig. 3b). In the next subsection, we will comment on possible correlations between these zero-temperature V vac results and the lattice-related contributions to the formation energy of V O at finite temperatures (i.e., the Gibbs free energy G * qh vac shown in Eq. (4)). The estimation of oxygen vacancy formation energies may depend strongly on the concentration of V O considered in the simulations due to the presence of short-and long-ranged interactions acting between the defects 2,8 . Figures 3a and 4 explicitly show this effect, as it is found that by decreasing the V O concentration the computed zero-temperature formation energy dramatically decreases for any arbitrary value of η . For instance, the estimated E vac for unstrained SrTiO 2.75 and SrTiO 2.94 amounts to 2.7 and 0.6 eV, respectively. This result suggests that short and middle-range interactions between oxygen vacancies are of repulsive type and thus the formation of V O clusters in STO thin films in principle is not likely to occur at low and moderate temperatures. Moreover, the E vac difference between equatorial and apical oxygen vacancies also depends critically on the concentration of defects. Specifically, according to our E vac results obtained for SrTiO 2.75 thin films in general it is more favourable to create apical V O than equatorial (Fig. 3a) whereas for SrTiO 2.94 thin films the tendency is just the opposite (Fig. 4b). The effect of epitaxial strain on E vac also varies as the concentration of V O changes. In particular, E vac increases both under compressive and tensile strains for SrTiO 2.94 thin films whereas for SrTiO 2. 75 it decreases under tensile strain.
How do these zero-temperature V O formation energy results compare with the available experimental data? In a recent paper, Rivadulla and collaborators have measured the enthalpy of oxygen vacancy formation for STO thin films as a function of epitaxial stress 34 . The authors have found that under both compressive and tensile strains the V O formation energy noticeably decreases. For instance, in the experiments the V O formation www.nature.com/scientificreports/ enthalpy decreases by ≈ 20 % ( ≈ 40 %) for a tensile (compressive) strain of 1% as compared to the unstrained case 34 . Therefore, the agreement between our zero-temperature E vac (η) results and the experimental observations is far from satisfactory. We note that this conclusion is independent of the V O concentration considered in our simulations, as shown by Figs. 3a and 4. In order to fundamentally understand the origins of such large discrepancies, and based on the fact that oxygen vacancies in oxide perovskites typically are created at high temperatures 7,25,34 , we proceeded to explicitly calculate V O formation free energies at finite temperatures (rather than for non-realistic T = 0 conditions). Before explaining our V O formation energy results obtained at T = 0 conditions, it is worth mentioning that in a recent work 33 another first-principles study on the V O formation energy of STO thin films has been reported. Zero-temperature E vac results analogous to ours are presented in 33 , however, the conclusions reported in that study are drastically different from the computational outcomes just described in this section. In particular, a systematic decrease in E vac has been predicted for either tensile or compressive strains, which is the opposite behaviour than what we have found here for SrTiO 2.94 thin films, for instance. Moreover, an intriguing correlation between the η-induced behaviour of E vac and the energy band gap of STO thin films ( E gap ) has been also suggested in work 33 . Based on our results enclosed in Figs. 2d and 4b, such a correlation is partially corroborated 21 . Nevertheless, in our calculations both quantities E vac and E gap increase, rather than decrease, under either tensile or compressive strains. We hypothesize that the likely reasons for such theoretical disagreements may be the neglection of characteristic STO structural motifs like polar and antiferrodistortive oxygen octahedral distortions in work 33   www.nature.com/scientificreports/ tion. Due to obvious computational limitations associated with the calculation of the full phonon spectrum of defective crystals, the dependence of G * qh vac on vacancy concentration could not be assessed. Figure 5 shows our G * qh vac results expressed as as function of temperature and epitaxial strain. Since here we are primarily interested in analyzing the joint effects of epitaxial strain and lattice thermal excitations on the formation energy of oxygen vacancies, the chemical potential entering Eq.(4) has been arbitrarily selected, without any loss of generality, to provide null G * qh vac values for the energy minimum determined at η = 0 and each temperature ("Free energy calculations" section).
We found that the η-dependence of the µ O -shifted formation energy of equatorial V O 's drastically changes as a result of considering T-induced lattice vibrations. For instance, at the highest analyzed temperature, T = 1000 K, G * qh vac (Eq) decreases by ≈ 30 % for a biaxial strain of −5 % and by ≈ 200 % for η = +5 % as compared to the T = 0 unstrained case (Fig. 5c). For an intermediate temperature of 500 K, the observed tendency is analogous to the one just described although the G * qh vac (Eq) differences with respect to the T = 0 unstrained case are slightly smaller (i.e., a reduction of ≈ 15 % and ≈ 150 % for η = −5 and +5 %, respectively- Fig.5b). Meanwhile, for apical oxygen vacancies G * qh vac (Ap) only changes moderately under tensile biaxial strains. The differences between the estimated G * qh vac as a function of T and η for apical and equatorial V O can be qualitatively understood in terms of the zero-temperature proxy V vac introduced in "Formation energy of oxygen vacancies at T = 0" section (Fig. 3b). In a recent theoretical paper 41 , it has been proposed that for positive V vac values, that is, V SrTiO 3−δ > V SrTiO 3 , lattice thermal excitations tend to facilitate the formation of oxygen vacancies. As it is observed in Fig. 3b, for equatorial vacancies V vac is positive and steadily increases under tensile biaxial strain; this outcome is agreeing with the large relative G * qh vac (Eq) decrease estimated for η = +5 % upon increasing temperature (Fig. 5). Meanwhile, for apical vacancies V vac is negative and may be small in absolute value under large tensile and compressive strains; this behaviour is consistent with the fact that under increasing temperature the corresponding relative G * qh vac (Ap) differences with respect to the T = 0 unstrained case only change slightly. Therefore, we corroborate the previously proposed qualitative correlation between the two quantities V vac and F qh vac ("Free energy calculations" section), which are computed at zero temperature and T = 0 conditions, respectively 41 .
How do these finite-temperature V O formation energy results compare with the experimental data reported in work 34 ? The answer is that although the agreement between theory and observations is not quantitative it can be regarded as qualitatively satisfactory. We recall that experimentally it has been determined that under both compressive and tensile biaxial strains oxygen vacancies can be created more easily. This behaviour is analogous to what we have predicted for equatorial V O 's, which in oxide perovskites correspond to the most representative class of anion positions (i.e., equatorial O sites are 50% more numerous than apical). Moreover, since the G * qh vac values estimated for equatorial V O 's under both tensile and compressive biaxial strains are smaller than those estimated for apical vacancies (by ≈ 30 and 20 meV, respectively), it is likely that to a certain extent vacancy www.nature.com/scientificreports/ ordering occurs in epitaxially strained STO thin films (as it has been experimentally shown for grain boundaries in bulk STO from scanning transmission electron microscopy measurements 75 ). On the down side, experiments indicate that it is more easy to create oxygen vacancies under compressive strain than under tensile strain 34 while our calculations predict the opposite trend (Fig. 5). Nonetheless, based on our computational E vac and G * qh vac results, it can be concluded that in order to reproduce the experimentally observed η-induced enhancement of V O formation with theoretical ab initio methods it is necessary to explicitly consider vibrational lattice thermal excitations in the calculations.
Zero-temperature activation energy for oxygen diffusion. The diffusion of O ions in oxide perovskites is a key parameter for the design of ionic-based devices 76 . In recent atomic force microscopy experiments performed by Iglesias et al. 37 , it has been shown that tensile biaxial strain produces a substantial increase in the diffusion of O ions in STO thin films. In particular, the room-temperature diffusion coefficient of oxygen atoms, D O , roughly increases by a factor of 4 upon a tensile biaxial strain of ≈ +2% 37 . For compressive tensile strains, on the other hand, the available experimental data is quite scarce. Nonetheless, measurements performed up to a η of ≈ −1 % appear to suggest an incipient reduction in D O 37 . On this regard, first-principles analysis of ionic transport properties may be very useful as calculations are free of the technical issues encountered in the experimental synthesis of epitaxially grown thin films and thus arbitrarily large tensile/compressive biaxial strains can be simulated. First-principles simulation of ionic diffusion processes, however, are neither exempt of some technical issues and shortcomings 59 . For instance, due to the intense computational expense associated with T = 0 simulations, most first-principles studies usually neglect temperature effects. In particular, zero-temperature calculations of ion-migration energy barriers typically are performed with the nudged elastic band (NEB) method 58 ("Zerotemperature activation energy for oxygen diffusion" section), in which (i) the initial and final diffusion positions of the vacancy and interstitial ions need to be guessed in the form of high-symmetry configurations rendering metastable states, and (ii) T-induced lattice excitations are totally neglected. Limitations of the NEB method for accurately determining ionic diffusion energy barriers and paths are well known and documented for some prototype fast-ion conductor materials (e.g., see works 59 Figure 6 shows our E a results obtained for (001) STO thin films by employing DFT NEB techniques ("Zerotemperature activation energy for oxygen diffusion" section). Two possible V O diffusion paths, namely, "Ap-Eq" (Fig. 6a) and "Eq-Eq" (Fig. 6c) where "Ap" and "Eq" stand for apical and equatorial O sites, have been considered in our simulations. In the former case, we obtain two different energy barriers, "Ap-Eq" and "Eq-Ap", due to the energy asymmetry between the two involved oxygen positions (Figs. 3 and 4). In consistent agreement with the available experimental data and previous DFT studies, we find that under tensile biaxial strain the energy barrier for V O diffusion is greatly reduced. For instance, at η ≈ +4 % we obtain that E a decreases with respect to the value estimated at zero strain (i.e., 0.55 eV) by ≈ 50 % and 15% for "Eq-Eq" (Fig. 6d) and "Eq-Ap" (Fig. 6b), respectively. (The V O diffusion energy barrier difference between cases "Eq-Ap" and "Ap-Eq" simply correspond to the zero-temperature V O formation energy difference between cases "Eq" and "Ap".) It is worth noting that our estimated zero-strain E a value of 0.55 eV is in very good agreement with the experimental V O diffusion energy barrier measured for bulk STO, E expt a ≈ 0.60 eV 78 . Upon compressive biaxial strain, we find that E a increases significantly and practically linearly with |η| (Fig. 6b,d). For instance, at η ≈ −4 % we predict that E a increases with respect to the zero-strain value of 0.55 eV by ≈ 45 % and 32% for "Eq-Eq" (Fig. 6d) and "Eq-Ap" (Fig. 6b), respectively. These results appear to be in agreement with the scarce experimental data that is available for compressive biaxial strains 37 but in apparent disagreement with previous DFT results reported by Al-Hamadany et al. 38 . The reasons for the disagreements between our theoretical NEB E a estimations and analogous computational results previously reported 38 could be the use of different exchange-correlation energy functionals and simulation supercells, which are known to affect the description of polar and antiferrodistortive oxygen octahedral distortions in functional oxide perovskites 47,79 . In order to fully test the reliability of our E a zero-temperature NEB results, we performed complementary ab initio molecular dynamics (AIMD) simulations in which lattice thermal excitations are fully taken into account and no particular V O diffusion path needs to be guessed 59 . Oxygen ionic diffusion at finite temperatures. Figure 7 encloses the MSD and D O results obtained from our T = 0 AIMD simulations for (001) STO thin films at η = ±3.6 % and zero strain ("Ab initio molecular dynamics simulations" section). For the η = 0 case, we estimate large diffusion coefficients of ∼ 10 −8 -10 −7 cm 2 s −1 at temperatures higher than 1000 K and a small V O diffusion energy barrier of 0.30 eV (Fig. 7b). The preexponential factor entering the corresponding D O Arrhenius formula ("Ab initio molecular dynamics simulations" section) amounts to 1.8 · 10 −6 cm 2 s −1 . It is important to note that the E a estimated by fully taking into account lattice thermal excitations is approximately 50% smaller than the corresponding value calculated with the NEB method by considering zero-temperature conditions. This computational outcome demonstrates the www.nature.com/scientificreports/ existence of an important interplay between lattice vibrations and V O diffusion, which in the case of STO thin films enormously facilitates ionic transport. It is also worth noting that the agreement between our zero-strain E a result obtained from AIMD simulations and the experimental diffusion energy barrier E expt a ≈ 0.60 eV 78 has considerably worsened as compared to the corresponding NEB estimation. Possible causes explaining such a disagreement could be the neglection of other types of crystalline defects that are relevant to ionic diffusion in our T = 0 calculations (e.g., dislocations 80 ), and the fact that the concentration of oxygen vacancies in our AIMD simulations ( ≈ 1.6 %) probably is higher than in the samples that were analyzed in the experiments.
For a tensile strain of +3.6 %, we find that the diffusion of oxygen ions is considerably enhanced as compared to the η = 0 case. In particular, we estimate high-T diffusion coefficients of ∼ 10 −7 cm 2 s −1 and a reduced O diffusion energy barrier of 0.17 eV (Fig. 7b). The value of the pre-exponential factor entering the corresponding D O Arrhenius formula ("Ab initio molecular dynamics simulations" section) is equal to 2.2 · 10 −6 cm 2 s −1 . The E a decrease induced by η = +3.6 % is about 50% of the zero-strain value, which is very similar to the relative variation estimated with NEB techniques for the same biaxial strain and "Eq-Eq" vacancy diffusion path ("Zerotemperature activation energy for oxygen diffusion" section). In this case, it is also concluded that the effects of lattice thermal excitations is to significantly enhance oxygen transport. As regards compressive biaxial strains, it is found that even at temperatures as high as 1500 and 2000 K the diffusion coefficient of oxygen atoms is nominally zero (Fig. 7a). This AIMD result is in qualitative agreement with the NEB calculations presented in the previous section, since in the latter case we found that E a increases almost linearly with |η| ("Zero-temperature activation energy for oxygen diffusion" section).  ) and (d). Labels "Eq" and "Ap" stand out for equatorial and apical oxygen vacancies, respectively. The colouring code for atoms in (a) and (c) coincides with that indicated in Fig. 1 www.nature.com/scientificreports/ Overall, the AIMD simulation results presented in this section confirm the correctness (at the qualitative level) of our NEB results reported in "Zero-temperature activation energy for oxygen diffusion" section, and demonstrate that lattice thermal vibrations have a significant enhancing effect on the diffusion of oxygen ions in (001) STO thin films. Interestingly, it is not always the case that lattice thermal excitations are found to promote ionic transport. For instance, in a recent systematic theoretical study on Li-based fast-ion conductors 59 the opposite trend has been demonstrated, namely, the energy barriers for ionic transport estimated from AIMD simulations in general are higher than those obtained with NEB methods. It is likely that the degree of anharmonicity of the non-diffusing lattice in the considered material, which determines the amplitude of the atomic fluctuations around the corresponding equilibrium positions, is directly related to the either enhancing or suppressing ionic diffusion effect mediated by the lattice excitations. Further quantitative investigations on this subject deserve future work.

Conclusions
We have presented a comprehensive ab initio study on the formation energy of oxygen vacancies and ionic diffusion properties of epitaxially strained (001) STO thin films, a class of functional materials with great fundamental and applied interests. The novelty of our work lies in the incorporation of lattice thermal excitations on the first-principles description of oxygen vacancies under varying epitaxial strain conditions. It has been demonstrated that in order to achieve an improved agreement with the experimental observations it is necessary to explicitly consider temperature-induced V O lattice effects in the theoretical calculations. By performing quasi-harmonic Gibbs free energy calculations, we have been able to qualitatively reproduce the nonmonotonic peak-like dependence of the V O formation enthalpy measured in STO thin films. Also, by performing ab initio molecular dynamics simulations we have been able to reproduce the qualitative η-driven oxygen ion diffusion trends observed in biaxially strained (001) STO samples. Generalization of our main conclusions to other technologically relevant oxide perovskite materials is likely although further experimental and computational works on the interplay between oxygen vacancies and epitaxial strain are necessary. We hope that the present study will stimulate new research efforts in this direction.

Data availability
The data that support the findings of this study are available from the corresponding author (C.C.) upon reasonable request. www.nature.com/scientificreports/ Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.