First principle electronic, structural, elastic, and optical properties of strontium titanate

We report self-consistent ab-initio electronic, structural, elastic, and optical properties of cubic SrTiO$_{3}$ perovskite. Our non-relativistic calculations employed a generalized gradient approximation (GGA) potential and the linear combination of atomic orbitals (LCAO) formalism. The distinctive feature of our computations stem from solving self-consistently the system of equations describing the GGA, using the Bagayoko-Zhao-Williams (BZW) method. Our results are in agreement with experimental ones where the later are available. In particular, our theoretical, indirect band gap of 3.24 eV, at the experimental lattice constant of 3.91 \AA{}, is in excellent agreement with experiment. Our predicted, equilibrium lattice constant is 3.92 \AA{}, with a corresponding indirect band gap of 3.21 eV and bulk modulus of 183 GPa.


I. INTRODUCTION
Strontium titanate (SrTiO 3 ) is one of the most studied oxides of the ABO 3 perovskite type structures, due to its great technological importance. Many interesting phenomena such as colossal magnetoresistance, high-T c superconductivity, multiferroicity, and ferroelectricity are observed in complex oxides. Since most of the interesting complex oxides have perovskite structure, SrTiO 3 is an ideal starting point for their study. It has been widely used for integration with other oxides into heterostructures. Those heterostructres show interesting properties such as thermoelectricity [1,2] and superconductivity [3,4]. Many new concepts of modern condensed matter and the physics of phase transitions have been developed while investigating this unique material [5][6][7]. SrTiO 3 has applications in the fields of ferroelectricity, optoelectronics and macroelectronics. It is used as a substrate for the epitaxial growth of high temperature superconductors. SrTiO 3 exhibits a very large dielectric constant. In comparison with SiO 2 , SrTiO 3 has almost two orders of magnitude higher dielectric constant and may as well offer a better replacement for SiO 2 in Si-based nanoelectronic devices (see Wilk et al. [8][and references therein]). SrTiO 3 has found usage in optical switches, grain-boundary barrier layer capacitors, catalytic activators, waveguides, laser frequency doubling, high capacity computer memory cells, oxygen gas sensors, semiconductivity, etc [9][10][11][12][13][14][15][16].
During the last few decades, the electronic, structural, elastic, and optical properties of SrTiO 3 (STO), as a model of ABO 3 perovskite, have been under intensive investigation both experimentally [17][18][19][20][21][22][23] and theoretically [9][10][11][12][13][24][25][26][27]. But, from a theoretical point of view, a proper description of its electronic properties is still an area of active research. Theoretical computations have had difficulty in predicting the correct band gap energy and other related electronic properties of SrTiO 3 from first principle. The density functional theory plus additional Couloumb interactions (DFT+U) formalism [28][29][30][31] has had good successes in obtaining correct energy bands and gaps of materials, but can only be applied to correlated and localized electrons, e.g., 3d or 4f in transition and rare-earth oxides. The hybrid functionals (for e.g., Heyd-Scuseria-Ernzerhof (HSE) hybrid functional [32][33][34]) has also been used in attempt to improve on the energy bands and band gaps of materials. This approach involves a range separation of the exchange energy into some fraction of nonlocal Hartree-Fock exchange potential and a fraction of local spin density approximation (LSDA) or generalized gradient approximation (GGA) exchange potential. We should note that this range separation is not universal. There is always a range separation parameter ω which varies between 0 and ∞. While it is reasonably clear that there exists a value of ω that gives the correct gap for a given system, this ω is not universal as it is always adjusted from one system to another [35,36]. For example, in HSE06 [33,37], ω = 0.11a −1 0 (a 0 is the Bohr radius) and in Perdew-Burke-Ernzerhof (PBEh) global hybrid [38], it is 25 % short-range exact exchange and 75 % shortrange PBE exchange. Even though the HSE functional, in most cases, accurately reproduces the optical gap in semiconductors, it severely underestimates the gap in insulators [36,39] and its band width in metallic systems is generally too large [36,[39][40][41]. The Engel and Vosko [42] (EV) GGA and the Tran and Blaha [43] modified Becke-Johnson (TB-mBJ) have also provided some improvements to the band gap of materials. For TB-mBJ, while the band gaps are considerably improved, the effective masses are severely underestimated [41]. In the case of the EV potential, the equilibrium lattice constants are far too large as compared to experiment and, as such, leads to an unsatisfactory total energy [44,45].
The theoretical underestimations of band gaps and other energy eigenvalues have been ascribed to the inadequacies of density functional potentials for the description of ground state electronic properties of semiconductors [18,19,24]. Also, other methods [46][47][48] that entirely go beyond the density functional theory (DFT) do not obtain the correct band gap values of most semiconductors without adjustment or fitting parameters [49,50]. This unsatisfactory situation is a key motivation for our work.
In this paper, we present a simple, yet robust, and abinitio method, based on self consistent solutions of the pertinent system of equations [51][52][53][54], that correctly predicts band gap values and related electronic properties of SrTiO 3 rigorously, from first principle, within the LCAO-GGA formalism. We also compute the structural, elastic and optical properties of SrTiO 3 .
The rest of this article is organized as follows. After this introduction in section I, the computational methods and details are given in section II. The results of our self-consistent calculations are presented and discussed in section III. We then summarize and conclude in section IV.

II. METHODS AND COMPUTATIONAL DETAILS
In the ground state, STO has the simple cubic (O 1 h -Pm3m) perovskite structure [55], with Sr atom sitting at the origin, Ti at the body center and three oxygen atoms at the three face centers [56] (see Fig. 1). We used the room temperature experimental lattice constant of 3.91 A [55,56].
Our ab initio, self consistent, nonrelativistic calculations employed a linear combination of atomic orbitals (LCAO) formalism and a generalized gradient approximation (GGA) potential. One may argue that relativistic effects are important for the description of SrTiO 3 . As was noted by Marques et al. [57], relativistic effects are only important for the description of the high-κ dielectric band structure. Their calculated, relativistic and nonrelativistic band structure for SrTiO 3 , up to an energy of 5 eV for the valence and conduction bands, respectively, are almost identical. Consequently, we expect a negligible relativistic correction for the band gap of SrTiO 3 .
Instead of assuming that a single trial basis set will yield the correct ground state charge density of the solid, the BZW method entails a minimum of three selfconsistent calculations with basis sets of different sizes and generally with different polarization functions, i.e., p, d, and f functions. The correct ground state is the one where all the occupied energies are at their minima. In practice, up to seven self-consistent calculations have been performed for some materials e.g, wurtzite ZnO [64]. The computations begin with a relatively small basis set that should not be smaller than the minimal basis set. The latter is defined as one that is just large enough to account for all the electrons in the atomic or ionic species present in the solid. The preliminary, self-consistent calculations of the electronic properties of the species pro-vide the wave functions that serve as input in the solid state calculations.
The first, self consistent calculation for the solid is performed with this small basis set (Calculation I) that is subsequently augmented with one orbital for the next self-consistent calculation (Calculation II). Depending on the angular symmetry of the added orbital, the size of the new basis set is larger than that of the previous one by 2, 6, 10, or 14 for s, p, d, and f functions, respectively. The occupied energies from calculations I and II are compared graphically and numerically. For the first two calculations, we found these occupied energies to be different for all the solids we have studied to date [65,66], including SrTiO 3 . The basis set for calculation II is then augmented in order to carry out self-consistent calculation III. Again, the occupied energies from calculations II and III are compared. This process of augmenting the basis set and of performing self-consistent calculations is continued until the occupied energies from a calculation, say N, are found to be identical to their corresponding ones from calculation (N+1), within computational uncertainties that are less than 50 meV. This perfect superposition of occupied energies from two consecutive calculations identifies the basis set for Calculation N as the optimal one, i.e, the smallest basis set that yields the lowest, occupied energies of the system. The attainement of this minima signifies that this basis set is verifiably complete for the description of the ground state of the system. Larger basis sets that include the optimal one do not lower any occupied energies from their values obtained with the optimal basis set.
As explained elsewhere [51,53], these larger basis sets do not lead to any changes in the ground state charge density or the Hamiltonian. However, larger basis sets that include the optimal one often lead to a lowering of some unoccupied energies. This lowering of some unoccupied energies cannot be ascribed to a physical interaction included in the Hamiltonian. Up to the optimal basis set, changes in the basis sets lead to changes in the charge density, the potential, and the Hamiltonian. Hence, changes in occupied and unoccupied energies, for self-consistent calculations with basis sets up to the optimal one, can be ascribed to a physical interaction embedded in the Kohn-Sham Hamiltonian. The system of equations in DFT totally determines changes in the occupied states. It also determines, at least in part, low-laying unoccupied states that are interacting with the occupied ones, up to the optimal basis set. For example, in wurtzite InN, the calculated dielectric functions agree with their corresponding experimental ones up to 5.5-6.0 eV [67]. Given that only direct transitions were taken into account in our dielectric functions calculations, this agreement indicates that the low-laying unoccupied bands were correctly determined. For larger basis sets that include the optimal one, the extra lowering of some unoccupied energies is a direct consequence of the Rayleigh theorem [53,58,59,65]. This theorem states that when an eigenvalue equation is solved with basis sets I and II, with set II larger than I and including I entirely, then the eigenvalues obtained with set II are lower than or equal to their corresponding ones obtained with basis set I [53,59].
The above process entails iterations for the equation giving the ground state charge density, with the iterations for the Kohn-Sham equation carried out for each choice of the basis set. Given that iterations for the Kohn-Sham (KS) equation involves the charge density (CD) equation, one could conclude that a single trial basis set calculation solves both equations self consistently. A problem with this view stems from the fact that any two such calculations, with different trial basis sets, lead to different, converged (i.e., self consistent) eigenvalues of the KS equation. The fundamental theorem of algebra suffices to guarantee that the two sets will be different if the basis sets have different numbers of basis functions utilized in the expansions. The question then arises as to which of the two sets of eigenvalues of the KS equation provides the physical description of the system under study. To answer such a question definitively and from first principle, the BZW method follows the process described above. Upon reaching the optimal basis set, not only the charge density, the potential, and the Hamiltonian no longer change (i.e., they have converged), but also the resulting, self-consistent eigenvalues of the KS equation have reached their respective minima, for the occupied states. In our understanding, to solve the system of equation self-consistently means obtaining converged eigenvalues (attainable with most arbitrary basis set) but also occupied eigenvalues that have reached their respective minima (attainable with BZW method).
In the above sense, the BZW method solves selfconsistently not just the Kohn-Sham equation, but also the equation giving the ground state charge density in terms of the wave functions of the occupied states. Further, in his Nobel lecture [68], Kohn noted the "density optimal" feature of the wave functions from correctly performed DFT calculations while those for the Hartree Fock approach are "total-energy optimal." Without a constrained search for the converged ground state, it is quite difficult to infer the basis set that yields the correct ground state charge density [69,70]. This point becomes clearer by noting that the reorganization of the cloud of valence electrons is drastically different for atomic or ionic species as compared to molecules or solids. For instance, single trial basis set calculations cannot make up for any deficiency in the angular symmetry of the functions, irrespective of the degree of convergence of the iterations of the Kohn-Sham equation. By correct ground state charge density, we mean the charge density that leads to the minima of all the occupied energies.
In this work, we utilize the electronic structure package from the Ames Laboratory of the US Department of Energy (DOE), Ames, Iowa [71]. We employ the generalized gradient approximation (GGA) potential given by Perdew and Wang [72]. We utilize sets of even tempered Gaussian functions with exponents from 0.12 to 10 5 to form the atomic wavefunctions. There are 15, 15, and 13 s, p, and d orbitals, respectively, for Sr, while 17, 17, and 15 s, p, and d orbitals, respectively, are used for both Ti and O. The charge fitting error using the Gaussian func-tions in the atomic calculation is about 10 −4 . Since the deep core states are fully occupied and are inactive chemically, the charge densities of the deep core states were kept the same as in the free atom. However, the core states of low binding energy were still allowed to fully relax, along with the valence states, in the self consistent calculations. The orbitals employed in the self-consistent calculations are between paranthesis for Sr (3d 4s 4p 4d 5s), Ti (3p 3d 4s) and O (2s 2p 3s), including some that are unoccupied in the free atoms or ions. These unoccupied orbitals are included in the self-consistent LCAO calculation to allow the restructuring of the electronic cloud, including possible polarization, in the crystal environment.
The Brillouin zone (BZ) integration for the charge density in the self consistent procedure is based on 56 special k points in the irreducible Brillouin zone (IBZ). The computational error for the valence charge is 5.3 x 10 −5 eV per valence electron. The self consistent potential converged to a difference of 10 −5 after several tenths of iterations. The energy eigenvalues and eigenfunctions are then solved at 161 special k points in the IBZ for the band structure. A total of 152 weighted k points, chosen along the high symmetry lines in the IBZ of SrTiO 3 , are used to solve for the energy eigenvalues from which the electron density of states (DOS) are calculated using the linear analytical tetrahedron method [73]. The partial density of states (pDOS) and the effective charge at each atomic sites are evaluated using the Mulliken charge analysis procedure [74]. We also calculated, the equilibrium lattice constant a o , the bulk modulus (B o ), the associated total energy and the electron and hole effective masses in different directions.
In calculating the lattice constant, we utilize a least square fit of our data to the Murnaghan's equation of state [75,76]. The lattice constant for the minimum total energy is the equilibrium one. The bulk modulus (B o ) is calculated at the equilibrium lattice constant.
The dielectric function ε(ω) = ε 1 (ω) + iε 2 (ω) can be calculated once the electronic wave functions and energies are known. The imaginary part of the dielectric function ε 2 (ω), from the direct interband transitions, is calculated using the Kubo-Greenwood formula [77]: where ω is the photon energy, P = −i ▽ is the momentum operator, Ω is the volume of the unit cell, ψ kn (r) and ψ kl (r) are the initial and final states, respectively, f ki is the Fermi distribution function for the i th states, and ǫ ki is the energy of the electron in the i th state. The real part of the dielectric function, ε 1 (ω), is obtained from the well-known Kramers-Kronig (KK) relation, where M indicates the principal value of the integral [78]. The real part of the optical conductivity, Re[σ(ω)], follows from above as

III. RESULTS AND DISCUSSION
The results of the electronic structure computations are given in Figs. 2 to 4. Figure 5 depicts the calculated total energy of STO, while Fig. 6 shows the optical spectra obtained using the optimal basis set from the electronic structure computations. Figures 1 and 4 have been drawn using the xcrysden [79]. We discuss the electronic structure, low laying conduction bands, and the effective mass in III A. The DOS is discussed elaborately in III B. The structural properties are presented in III C, while III D deals with the calculated, optical properties.

A. The Electronic Structure, Band Gap, Low-energy Conduction Band and Effective Mass
The electronic structure of the valence and the low energy conduction states determine the band gap and other important properties of materials. Table I shows that our ab initio, first principle method yielded an indirect band gap of 3.24 eV at the experimental lattice constant of 3.91Å (see Fig. 2) and an indirect band gap of 3.21 eV at the calculated equilibrium lattice constant of 3.92Å. This table contains several, previous results from calculations using LDA, GGA or hybrid potentials. Table II contains the calculated energies at some high symmetry points in the BZ. These energies are provided for future comparison with experimental and theoretical findings. Our calculated band structure (see Fig. 2) resembles that of the parent TiO 2 system. Fig. 2 also shows that the top of the valence band is at the L point.
The effective mass is one of the main factors determining the transport properties, the Seebeck coefficient, and electrical conductivity of materials. The calculated electron effective masses at the bottom of the conduction band along the Γ -L, Γ -X, and Γ -K directions, respectively, are 0. 68   vations of Mattheiss and co-workers [80,81]. Our calculated effective masses are in excellent agreement with the detailed effective mass values as reported by Mattheiss [80][and references therein] and the relativistic computational results of Marques et al. [57].
There is a significant O 2p -Ti 3d hybridization in the valence bands. As shown in Fig. 2 results of Jiangni et al. [12]. They are formed by the hybridization between O 2p and Ti 3d, with very little contribution from the Ti 3p and Sr 4p (see Fig. 3(b)). Two of the bands at the Γ point are triply degenerate: 1 Γ 15 (-2.75 eV) and 2 Γ 15 (-0.25 eV), while the third band Γ 25 (-1.23 eV) is non-degenerate (see Table II  Ref. [10]. b Ref. [6]. c Ref. [12]. d Ref. [57]. e Ref. [82]. f Ref. [83]. g Ref. [84]. h Ref. [49]. i Ref. [9]. j Ref. [5]. k Ref. [13]. l Ref. [56]. m Ref. [85]. n Ref. [86]. position of the O 2s to be about 17 eV below E F . Our calculated value of -17.85 eV in Fig. 3(a) is close to the experimental one. In particular, our result does not underestimate this core state position as is often the case in GGA calculations. The conduction bands (CB), immediately above the Fermi level (low energy conduction band), are dominated by threefold degenerate Ti 3d t 2g orbitals which hybridize with O 2p and O 2s. The two-fold degenerate Ti 3d e g states have some hybridization from all other orbitals except Ti (4s and 3p), Sr (4p and 3d). The energy eigenvalue in the lowest conduction bands, at the X point, are practically the same as that at the Γ point, resulting in the observed, minimal dispersion in the conduction band between Γ and X points. This feature is apparent from our energies at the high symmetry points in Table II. At the Γ point, energies associated with the lowest-laying conduction bands are: 1 Γ 25 ′ (3.24 eV), 1 Γ 12 (4.84 eV), 2 Γ 12 (6.10 eV), 3 Γ 1 (6.75 eV) and 2 Γ 25 ′ (8.92 eV). The calculated, low energy conduction bands in Fig. 2 are quite different from those of previous studies. Figure 2 shows that the lowest conduction bands are degenerate at the Γ and X points along the [100] and equivalent directions. The electronic structure in Fig. 2 was calculated using the experimental lattice constant. We further examined whether or not the position of the shallow minimum in the lowest conduction band depends on the value of the lattice constant by using several values of the lattice constant around the experimental one. Even though, the band gap value changed from 3.26 to 3.17 eV, there was no appreciable change in the depth of the shallow minimum in the lowest conduction band. We recall that the gap is 3.21 eV for our calculated, equilibrium lattice parameter.   [90]. c Ref. [11]. d Ref. [91]. e Ref. [85]. f Ref. [86]. g Ref. [55,56,92,93]. h Ref. [56,92,93].

B. Densities of States, Electron Distribution and Chemical Bonding
Figures 3(a) and 3(b) exhibit the total (DOS) and related partial (pDOS) densities of states, respectively. Figure 4 shows the contour plot of the distribution of the electron charge density of SrTiO 3 . As can be seen from Fig. 4 [94,95] with corresponding charge densities of 0.2 and 0.67 -0.90 e/Å 3 , in that order [94][95][96][97]. The charge density distribution around the O atom with respect to the horizontal Ti -O -Ti line is elongated in the direction along the Ti -O covalent bond in agreement with room temperature experimental charge density distribution of SrTiO 3 reported by Ikeda et al. [96]. This anisotropic charge density distribution is ascribed to the rotational mode of the Ti -O 6 octahedron by experiment [94,95,98].
From our calculated DOS (see Fig. 3(a)), it can be inferred that the onset of absorption is quite sharp and it starts at about 3.24 eV. It exhibits a fine structure at 3.6 eV, with a shoulder at 4.50 eV. This picture is consistent with the experimental results of Cohen and Blunt [99] and the theoretical findings of Perkins and Winter [100] of a relatively sharp absorption edge in the optical measurement of SrTiO 3 . In the calculated DOS of the low laying conduction bands, sharp peaks appear at 4.70 eV, 5.72 eV, and 7.05 eV. Relatively broad peaks are found at 8.33 eV, 10.97 eV and 12.75 eV. Our computed peaks are in basic agreement with experimental findings of Cardona [101] and Braun et al. [102]. For the valence bands DOS, we calculated peaks at -0.20 eV, -0.72 eV, -1.14 eV, -1.83 eV, -2.86 eV, -4.06 eV, -4.50 eV, -14.30 eV, -4.61 eV, and -17.53 eV. The peaks in the valence band DOS are all sharp. Our calculated electronic structure is in agreement with scanning transmission electron microscopy, vacuum ultraviolet spectroscopy, and spectroscopic ellipsometry measurement of Van Benthem et al. [103].
Our calculated band gap value of 3.24 eV, from L to Γ, is practical the same as the experimental one. In general, other theoretical calculations obtained values that are up to 1.1 eV smaller. The source of the small band gap values was believed to be due to the pushing up of the top of the valence band dominated by Ti 3d and O 2p states to higher energies [22]. According to our findings, it rather appears to be the extra lowering of the conduction bands that produces GGA (or LDA) band gaps that are more than 1.1 eV smaller than the experimental values, if LCAO type computations do not search for and utilize an optimal basis set. Such a basis set is verifiably converged for the description of occupied states [53,60,61]. We recall that Kohn and Sham [104], in their original paper, explicitly stated the need to solve self consistently the pertinent system of equations defining LDA. The BZW method, as explained above, rigorously solves the system of equations in the sense explained above. Single trial basis set calculations also involve both the KS and charge density equations. The major difference resides in the fact that these calculations do not generally entail changes in the basis functions beyond those of the expansion coefficients.

C. Structural Properties
The total energy versus the lattice constant data are shown in Fig. 5. The data fit well to the Murnaghan equation of state (EOS). The calculated equilibrium lattice constant is 3.92Å and the bulk modulus, B o , is 183.45 GPa.
The experimently reported lattice constants are in the range 3.89 to 3.92Å [55,56,86,92] and the bulk modulus lays in the range 174 to 183 GPa [56,86,92]. In Table III, we show our calculated equilibrium lattice constant and bulk modulus in comparison with experimental and other theoretical results. Both our calculated lattice constant and bulk modulus agree well with corresponding, experimental ones, respectively.

D. Optical Properties
The plot of the dispersive (ε 1 (ω)) and absorptive (ε 2 (ω)) dielectric functions are shown in Figs. 6(a) and 6(b), respectively, while the optical conductivity (σ(ω)) profile is in Fig. 6(c). All reported spectra have been calculated without any broadening and, may have more features than experimental ones. Our calculated, dielectric spectra are in good agreement with the experimental measurements [93,101]. The calculated optical spectra only included the direct inter band transitions. The fundamental absorption edge E o , which is also a measure of the optical gap, was found to be 3.24 eV from the calculation, as per the insert of Fig. 6(b). Our computed direct gap of 3.43 eV is in agreement with the experimental value of 3.40 eV [93,101]. Our calculated ε 1 (ω = 0) at zero frequency equals 4.75 (cf. Fig. 6(a)). It compares well with the experimental value of 4.92 measured by Braun et al. [102]. Our dielectric spectra resemble that of BaTiO 3 of Bagayoko et al. [105]. This observation also holds for the data of Cardona [101] and Baurerle et al. [106]. Both the experimental and our calculated results show that the direct optical gap is larger than the smallest indirect band gap. Figure 6(c) shows the optical conductivity σ(ω) of SrTiO 3 . As per the insert, it also shows that the fundamental absorption edge starts at 3.24 eV. The positions of the peaks (without any rigid shift) are in agreement with experimental findings.

IV. SUMMARY AND CONCLUSION
We have performed first principle, ab initio calculations of the electronic, structural, elastic, and optical properties of bulk SrTiO 3 in the cubic phase using GGA potential and the BZW method. Our calculated results, without any adjustment or corrections, show good agreement with experimental data.
The agreement of our calculated band gaps (3.21 and 3.24 eV) and electron effective masses with corresponding, experimental values is significant. Some calculations with adjustable parameters can lead to the correct band gap; but they generally do not yield the correct curvature of the conduction band around its minimum-as given by the electron effective masses. Similarly, the agreement between the peaks in the calculated density of states with corresponding, experimental ones denote the correct description of the relative location of the bands. This result is confirmed by our reproduction of the measured features of the dielectric functions, the imaginary part of which was obtained using only direct transitions between occupied and unoccupied bands. Our calculated equilibrium lattice constant (3.92Å) and bulk modulus (183.45 GPa) also agree with corresponding, experimental findings.