text:
# Nature of the insulating phases in the half-filled ionic Hubbard model

## I Introduction

For more than two decades the correlation induced metal-insulator transition and its characteristics has been one of the challenging problems in condensed matter physics . This metal-insulator transition is often accompanied by a symmetry breaking and the development of long range order . In one dimension this ordering can only be related to the breaking of a discrete symmetry. Examples include commensurate charge density waves (CDWs) and Peierls dimerization (bond-order wave, BOW) phenomena. In contrast, the transition into the Mott insulating (MI) phase in one dimension is not connected with the breaking of a discrete symmetry . In the MI phase the gapped charge degrees of freedom are uniformly distributed in the system, while the gapless spin degrees of freedom are described by an effective $`S=1/2`$ Heisenberg chain .

Due to the different symmetries of the CDW, BOW, and MI phases it is natural to expect that these phases are mutually exclusive. The extended Hubbard model at half-filling with an on-site ($`U`$) and a nearest neighbor ($`V`$) Coulomb repulsion provides a prominent example with a transition from a MI to a CDW insulator in the vicinity of the $`U=2V`$ line in the phase diagram . Remarkably, the transition at weak coupling may involve an intermediate BOW phase . A similar phase-diagram structure was recently also discovered for the Holstein Hubbard model . The tendency towards BOW order is even more profound for the Hubbard model with an explicit bond-charge coupling where the CDW and MI phases are often separated by a long range ordered Peierls dimerized phase .

In recent years particular attention has been given to another example for an extension of the Hubbard model which includes a staggered potential term . The corresponding Hamiltonian has been named the “ionic Hubbard model” (IHM); in one dimension it is given by

$`H=`$ $``$ $`t{\displaystyle \underset{i,\sigma }{}}(1+(1)^i\delta )(c_{i\sigma }^{}c_{i+1\sigma }^{}+H.c.)`$ (1)
$`+`$ $`U{\displaystyle \underset{i}{}}n_i^{}n_i^{}+{\displaystyle \frac{\mathrm{\Delta }}{2}}{\displaystyle \underset{i,\sigma }{}}(1)^in_{i\sigma }^{},`$ (2)
where $`c_{i\sigma }^{}`$ creates an electron on site $`i`$ with spin $`\sigma `$ and $`n_{i\sigma }^{}=c_{i\sigma }^{}c_{i\sigma }^{}`$. $`\mathrm{\Delta }`$ is the potential energy difference between neighboring sites, and $`\delta `$ a Peierls modulation of the hopping amplitude $`t`$. In the limit $`\mathrm{\Delta }=\delta =0`$, Eq. (2) reduces to the ordinary Hubbard model, the limit $`\mathrm{\Delta }=0`$ and $`\delta >0`$ is called the Peierls-Hubbard model, and the limit $`\mathrm{\Delta }>0`$ and $`\delta =0`$ is usually referred to as the IHM. In the following, we will focus mainly on the effect of the on-site modulation $`\mathrm{\Delta }`$, so we implicitly assume $`\delta =0`$ except where stated otherwise.

The IHM was first proposed and discussed almost 20 years ago in the context of organic mixed-stack charge-transfer crystals with alternating donor ($`D`$) and acceptor ($`A`$) molecules ($`\mathrm{}D^{+\rho }A^\rho \mathrm{}`$) . These stacks form quasi-1D insulating chains, and at room temperature and ambient pressure are either mostly ionic ($`\rho 1`$) or mostly neutral ($`\rho 0`$). However, several systems undergo a reversible neutral to ionic phase transition i.e. a discontinuous jump in the ionicity $`\rho `$ upon changing temperature or pressure . Later the IHM has been used in a similar context to describe the ferroelectric transition in perovskite materials such as BaTiO3  or KNbO3 .

The very presence of at least on transition in the ground state phase diagram of the half-filled IHM model is easily traced by starting from the atomic limit . For $`t=0`$, it is obvious that at $`U<\mathrm{\Delta }`$ the ground state of the IHM has two electrons on the odd sites, and no electrons on the even sites corresponding to CDW order with maximum amplitude. On the other hand, for $`U>\mathrm{\Delta }`$ each site is occupied by one electron with infinite spin degeneracy. Thus, for $`t=0`$ a transition occurs at a critical value $`U_c=\mathrm{\Delta }`$. This transition is expected to persist for finite hopping amplitudes $`t>0`$.

A renewal of interest in the IHM started with the bosonization analysis of Fabrizio et al. (FGN) , where a two-transition scenario for the ground-state phase diagram of the 1D IHM was proposed. The key features of the FGN theory are the presence of an Ising-type transition from a CDW band-insulator phase at $`U<U_{ch}^c`$ into a BOW phase at $`U_{ch}^c<U<U_{sp}^c`$, and a continuous Kosterlitz-Thouless like transition into a MI phase at $`U>U_{sp}^c`$. In this scenario the charge gap vanishes only at $`U=U_{ch}^c`$ and the system might be “metallic” at this point. The second transition at $`U=U_{sp}^c`$ is connected with the closing of the spin gap, which is finite for all $`U<U_{sp}^c`$ and vanishes for $`U>U_{sp}^c`$. Thus, the bosonization phase diagram essentially supports the ”exclusion principle” of the ground states.

Later on various attempts based on numerical tools have been performed to verify the FGN phase diagram for the half-filled IHM. In particular, exact diagonalization , valence bond techniques , quantum Monte Carlo , and DMRG  was used. Unfortunately, conflicting results have so far been reported in these studies regarding the nature of the transition and the insulating phases, the possibility of two rather than one critical point, or the appearance of BOW order.

Given the numerous unresolved issues we reinvestigate the ground-state properties of the IHM using the exact diagonalization Lanczos technique and the DMRG method. We verify the presence of at least one transition at a critical coupling $`U_c(\mathrm{\Delta })`$ from a band-insulator (BI) to a correlated insulator (CI) phase. On finite systems the transition originates from a ground-state level crossing with a change of the site-parity eigenvalue, which implies the vanishing of the optical excitation gap at $`U_c`$; Our DMRG results show that the spin and charge gaps remain nevertheless finite and equal at the transition. Above $`U_c`$ the charge and spin gaps split, the charge gap increases, while the spin gap decreases and we identify long-range BOW order with a spontaneous site-inversion symmetry breaking. The existence of a second transition is not unambiguously resolved within the accuracy of our DMRG data. Yet, the scaling of the BOW order parameter changes qualitatively with increasing $`U`$ indicative for a possible second smooth transition point where the spin gap closes.

We show that at $`U<U_c(\mathrm{\Delta })`$ the CDW-band insulator phase is realized. In this phase BOW and spin density wave (SDW) correlations are strongly suppressed, and the spin and charge gaps are equal and finite. The characteristic feature of the ground-state phases for $`U>U_c(\mathrm{\Delta })`$ is the coexistence of long range CDW order with either a long range BOW or algebraically decaying BOW and SDW correlations.

## II Symmetry analysis

A good starting point for understanding the existence of a phase transition in the IHM is to study the symmetry of the model manifestly seen in the limiting cases $`U\mathrm{\Delta },t`$ and $`U\mathrm{\Delta },t`$. The IHM is invariant with respect to inversion at a site and translation by two lattice sites. If we denote the site inversion operator by $`P`$, defined through

$$Pc_{i\sigma }^{}P^{}=c_{Li\sigma }^{}\mathrm{for}i=0,\mathrm{},L1,$$
(3)
and $`\widehat{T}_j`$ for a translation by $`j`$ sites, then any nondegenerate eigenstate $`|\psi _n`$ of $`H`$ must obey $`P|\psi _n=\pm |\psi _n`$ and $`\widehat{T}_2|\psi _n=|\psi _n`$. Because $`[H,\widehat{T}_1]0`$, any non-degenerate eigenstate $`|\psi _n`$ of $`H`$ is not an eigenstate of $`\widehat{T}_1`$.

For the half-filled Hubbard model ($`\mathrm{\Delta }=\delta =0`$) the ground state has $`P=+1`$ only for $`U=0`$, and $`P=1`$ for any $`U>0`$ . However, in the IHM the phase transition from a renormalized BI to a CI occurs at some finite $`U_c>0`$. This suggests that the parity of the ground state remains even not only for $`U=0`$, but for all $`U<U_c`$. At $`U_c`$, a ground-state level crossing occurs on finite chains, as confirmed by exact diagonalization studies (see below), connected with a site-parity change.

For $`U=0`$ the ground state at half-filling is a CDW-BI. The alternating potential defines two sublattices, doubling the unit cell and opening up a band gap $`\mathrm{\Delta }`$ for $`U=0`$ at $`k=\pm \pi /2`$. The elementary spectrum consists of particle-hole excitations over the band gap. The charge ($`\mathrm{\Delta }_C`$) and spin ($`\mathrm{\Delta }_S`$) excitation gaps are equal: $`\mathrm{\Delta }_C=\mathrm{\Delta }_S=\mathrm{\Delta }`$. We consider a system to be a BI when the criterion $`\mathrm{\Delta }_S=\mathrm{\Delta }_C`$ holds, where the spin and the charge gaps are given by

$`\mathrm{\Delta }_S`$ $`=`$ $`E_0(N=L,S_z=1)`$ (4)
$``$ $`E_0(N=L,S_z=0),`$ (5)
$`\mathrm{\Delta }_C`$ $`=`$ $`E_0(N=L+1,S_z=1/2)`$ (6)
$`+`$ $`E_0(N=L1,S_z=1/2)`$ (7)
$``$ $`2E_0(N=L,S_z=0),`$ (8)
respectively. $`E_0(N,S_z)`$ is the ground-state energy, $`L`$ the system length, $`N`$ the number of electrons, and $`S_z`$ the $`z`$-component of the total spin. As we show below the BI phase is realized in the ground state of the IHM at $`U<U_c`$.

In the strong-coupling limit $`U\mathrm{\Delta },t`$, the low-energy physics of the IHM is described by the following effective Heisenberg spin model

$$H_{eff}=J\underset{i}{}𝐒_i𝐒_{i+1}+J^{}\underset{i}{}𝐒_i𝐒_{i+2}.$$
(9)
In Eq. (9) the exchange couplings are given by

$`J`$ $`=`$ $`{\displaystyle \frac{4t^2}{U}}\left[{\displaystyle \frac{1}{1x^2}}{\displaystyle \frac{4t^2}{U^2}}{\displaystyle \frac{1+4x^2x^4)}{(1x^2)^3}}\right],`$ (10)
$`J^{}`$ $`=`$ $`{\displaystyle \frac{4t^4}{U^3}}{\displaystyle \frac{(1+4x^2x^4)}{(1x^2)^3}},`$ (11)
where $`x=\mathrm{\Delta }/U`$. This result (9) implies that in the strong-coupling limit of the IHM the low-energy physics is qualitatively similar to that of the Hubbard model, with modified exchange coupling constants $`J`$ and $`J^{}`$. For next-nearest neighbor couplings $`J^{}<0.24J`$ the spin gap vanishes . The coupling constants (11) satisfy the condition $`J^{}<0.24J`$ at least for $`U>3.6t`$ for $`\mathrm{\Delta }t`$ and $`U>3.6\mathrm{\Delta }`$ for $`\mathrm{\Delta }>t`$.

The effective spin model (9) is invariant with respect to translations by one lattice spacing, whereas the original IHM is invariant only with respect to translations by two lattice spacings. However, the doubling of the unit cell is ensured due to the charge degrees of freedom. For the standard Hubbard model at arbitrary $`U0`$ the number of doubly occupied sites $`D`$ in the ground state is finite. The exact Bethe-ansatz solution tells that $`D`$ scales as $`(t/U)^2`$ in the strong coupling limit and is given by

$$D=\underset{i}{}n_i^{}n_i^{}NA(t/U)^2[1+𝒪\left((t/U)^2\right)],$$
(12)
where $`A=4\mathrm{ln}2`$. Contrary to the Hubbard model, where doublons are equally distributed on all sites of the system, the non-equivalence of sites in the IHM leads to different probabilities for finding a doublon on even or odd sites. Since doublons are spin singlets, their distribution is not influenced by the spin fluctuations. Since the energy of doublons on even and odd sites differ in $`\mathrm{\Delta }`$ and assuming the scaling for the density of doublons as in Eq. (12), one easily obtains for the amplitude of the ionicity induced CDW in the strong coupling limit

$`{\displaystyle \frac{1}{N}}\left(D_{\mathrm{𝑜𝑑𝑑}}D_{\mathrm{𝑒𝑣𝑒𝑛}}\right)A_1{\displaystyle \frac{t^2}{U^2}}\left[{\displaystyle \frac{1}{(1x)^2}}{\displaystyle \frac{1}{(1+x)^2}}\right]`$ (13)
$`=4A_1{\displaystyle \frac{t^2}{U^2}}{\displaystyle \frac{x}{(1x^2)^2}}[1+𝒪\left((t/U)^2\right)]`$ (14)
where $`A_1`$ is a constant of order unity. Thus, although the effective spin Hamiltonian has a higher symmetry than the original model from which it was derived, the translational symmetry of the IHM is recovered due to the long range CDW pattern arising from the staggered doublon and holon distribution.

## III Exact diagonalization results

In order to explore the nature of the spectrum and the phase transition, we have diagonalized numerically small systems by the Lanczos method  similarly to earlier exact diagonalization calculations . The energies of the few lowest eigenstates were obtained for finite chains with $`L=4n`$ and periodic boundary conditions or $`L=4n+2`$ with antiperiodic boundary conditions.

We first analyze short chains; for chain lengths $`L16`$ finite-size effects do not change the qualitative behavior discussed below. In Fig. 1, the lowest eigenenergies of the IHM for $`\mathrm{\Delta }=0.5t`$, $`L=8`$ and periodic boundary condition are shown as a function of $`U`$. At $`U=1.3t`$, a level crossing of the two lowest eigenstates occurs. A non-degenerate eigenstate of the IHM has a well defined site-parity, so a ground-state level-crossing transition necessarily corresponds to a change of the site-parity eigenvalue.

For $`U=0`$, the IHM is easily diagonalized in momentum space by introducing fermionic creation operators $`\gamma _{k\sigma b}^{}`$ with a band index $`b=1,2`$ for the lower and upper bands, respectively, with the dispersion $`E_{1/2}(k)=\pm \sqrt{4\mathrm{cos}^2(k)+(\mathrm{\Delta }/4)^2}`$ for momenta $`\pi /2<k\pi /2`$. For $`U=0`$ the first two degenerate excited states at half-filling always have negative site parity, because the ground state has $`P=+1`$, and the operator $`\gamma _{q\sigma 2}^{}\gamma _{q\sigma 1}^{}`$ with $`q=\pi /2`$ obeys

$$P\gamma _{q\sigma 2}^{}\gamma _{q\sigma 1}^{}=\gamma _{q\sigma 2}^{}\gamma _{q\sigma 1}^{}P.$$
(15)
The first two excited states shown in Fig. 1 are the spin singlet ($`S=0`$, $`S_z=0`$) and triplet excitations ($`S=1`$, $`S_z=0`$), created from the ground state by applying the operators

$`{\displaystyle \frac{1}{\sqrt{2}}}\left(\gamma _{q2}^{}\gamma _{q1}^{}\gamma _{q2}^{}\gamma _{q1}^{}\right),`$ (16)
$`{\displaystyle \frac{1}{\sqrt{2}}}\left(\gamma _{q2}^{}\gamma _{q1}^{}+\gamma _{q2}^{}\gamma _{q1}^{}\right),`$ (17)
respectively. Thus both excited states have total momentum $`k_{tot}=0`$ and negative site parity. For $`U>0`$, these degenerate excited states split in energy. Exact diagonalization of finite IHM rings therefore identifies one critical $`U_c>0`$, separating a BI with $`P=+1`$ at $`U<U_c`$ from a CI with $`P=1`$ for $`U>U_c`$.

## IV DMRG results

### A Excitation gaps

In order to access the transition scenario in the long chain-length limit, we have studied chains up to L=512 using the DMRG method . The fact that the transition at $`U_c`$ is connected to a change in inversion symmetry requires some caution when open boundary conditions (OBC) are used in DMRG studies. For OBC and $`L=2n`$ the IHM is not reflection symmetric at any site. Thus, the ground state does not have a well defined site parity, and the level-crossing transition is absent. To overcome this problem, one might try to use chains with OBC and an odd number of sites $`L=2n+1`$, since the Hamiltonian in this case is reflection symmetric with respect to the site $`i_c`$ in the center of the chain, and a site inversion operator is well defined by

$$Pc_{i_c\sigma }^{}P^{}=c_{L1i_c\sigma }^{}i=0,\mathrm{},L1.$$
(18)
To test whether this is an improved choice we have calculated the site parity of the ground state for $`U=0`$ analytically for different chain lengths $`L=2n+1`$ and found

$$P|\psi _0=(1)^n|\psi _0.$$
(19)
On the other hand, if one extends the idea of Gidopoulos et al.  for the determination of the site parity to chains with $`L=2n+1`$ for $`Ut`$, one obtains

$$P|\psi _0=(1)^{\left[_{m=1}^{L1}m\right]}|\psi _0=(1)^n|\psi _0.$$
(20)
Thus, the parity eigenvalue of the ground state is the same at $`U=0`$ and $`Ut`$ for a given chain length, and no level crossing occurs. Due to the fact that the sharp transition at a well defined $`U_c`$ does not exist in the finite-chain results for OBC, the extrapolation is a rather subtle problem, since a sharp transition feature has to be identified from the extrapolation of smooth curves. This requires the use of quite long chains in the critical region.

In Fig. 2 extrapolated results are shown for the spin and charge gaps $`\mathrm{\Delta }_S`$ and $`\mathrm{\Delta }_C`$, respectively. Calculations were performed with OBC for chains of lengths $`L=\{30,40,50,60\}`$, and additionally up to $`L=512`$ in the transition region around the estimated $`U_c`$. We assume a scaling behavior of $`\mathrm{\Delta }_C`$ and $`\mathrm{\Delta }_S`$ of the form

$$\mathrm{\Delta }_i(L)=\mathrm{\Delta }_i^{\mathrm{}}+\frac{A_i}{L}+\frac{B_i}{L^2},$$
(21)
where $`i\{S,C\}`$. The extrapolation for $`L\mathrm{}`$ is then performed by fitting this polynomial in $`1/L`$ to the calculated finite-chain results. We note that different finite-size scaling formulas were proposed in the literature mainly when periodic or antiperiodic boundary conditions were used .

As can be seen from the main plot in Fig. 2, extrapolating the results for $`L=\{30,40,50,60\}`$ does indeed not give a sharp transition behavior. As illustrated in the inset, adding results for $`L`$ up to $`512`$ in the critical region changes the picture considerably. Within numerical accuracy the charge and spin gaps remain equal and finite up to a critical $`U_c2.1t`$, where a sharp kink for $`\mathrm{\Delta }_C`$ is observed. Importantly, $`\mathrm{\Delta }_C`$ does not close at the critical point. We emphasize that the magnitude of $`\mathrm{\Delta }_C`$ and $`\mathrm{\Delta }_S`$ at the transition point is sufficiently larger than our numerical uncertainty in the finite size scaling analysis and therefore allows for a safe conclusion. $`\mathrm{\Delta }_C=\mathrm{\Delta }_S>0`$ at the transition is in fact not in conflict with an underlying ground-state level crossing. If the ground states of the different site-parity sectors become degenerate, the only rigorous consequence is the closing of the optical excitation gap. The selection rules for optical excitations allow only for transitions between states of different site-parity. Furthermore, optical transitions occur within the same particle number sector. The optical gap is therefore by definition distinct from the charge gap Eq. (2) which involves the removal or the addition of a particle. The critical point $`U_c`$ of the IHM has the remarkable peculiarity that the optical gap closes while $`\mathrm{\Delta }_C`$ remains finite. Above $`U_c`$ the charge and spin gaps split indicating that the corresponding insulating phase is no longer a BI. $`\mathrm{\Delta }_S`$ continuously decreases with increasing $`U`$ and becomes unresolvably small above $`U2.5t`$ within the achievable numerical accuracy.

The result, that $`\mathrm{\Delta }_C`$ and $`\mathrm{\Delta }_S`$ remain finite and equal at the transition is in agreement with the data obtained by Qin et al. . These authors performed DMRG calculations for the IHM with $`\mathrm{\Delta }=0.6t`$ for open chains up to $`L=600`$ sites. They observed a surprising non-monotonic scaling behavior of $`\mathrm{\Delta }_S`$ with $`L`$ for values of $`U`$ close to the critical $`U_c`$, i.e. for chain lengths $`L>300`$ $`\mathrm{\Delta }_S`$ started to increase again. It remains unclear whether this is due to loss of DMRG accuracy with increasing chain lengths and keeping a fixed number of states in the DMRG algorithm. In contrast, our data always show a monotonous scaling with $`1/L`$. DMRG calculations for the IHM with $`\mathrm{\Delta }=t`$ have also been performed by Takada and Kido for chains up to $`L=400`$ sites . The authors interpret their results in the region close to $`U_c`$ in favour of a two-transition scenario similar to that of FGN . However, as we will show below the accuracy of the currently available DMRG data is not enough to provide a stringent argument in favour of this interpretation.

For comparison we show in Fig. 3 the spin and charge gaps versus $`U`$ in the Peierls Hubbard model. As we observe this model is distinctly different from the IHM. This is also a band-insulator at $`U=0`$, but in contrast to the IHM has $`\mathrm{\Delta }_C>\mathrm{\Delta }_S>0`$ for any value $`U>0`$, i.e. the phase transition from the Peierls band-insulator to the correlated insulator occurs at $`U_c=0`$. So although the Peierls and the ionic BI for $`U=0`$ similarly possess an excitation gap at the Brillouin zone boundary, applying a Coulomb $`U`$ leads to distinctly different behavior in both cases. The origin of the different behaviors must be traced to the fact that the Hubbard interaction and the ionic potential compete locally on each site, while the Peierls modulation of the hopping amplitude tends to move charge to the bonds between sites, thereby avoiding conflict with the Hubbard term.

The spin-Peierls physics of the Peierls Hubbard model at large $`U`$ evolves smoothly with decreasing $`U`$ into the physics of a spin-gapped CI-BOW state in the weak-coupling limit, which is characterized by long-range staggered bond-density correlations

$`g_b(r)`$ $`=`$ $`{\displaystyle \frac{1}{L}}{\displaystyle \underset{i}{}}\psi _0|b(i)b(i+r)|\psi _0,`$ (22)
$`b(i)`$ $`=`$ $`{\displaystyle \underset{\sigma }{}}(c_{i\sigma }^{}c_{i+1\sigma }^{}+H.c.).`$ (23)
For $`Ut`$ the low-energy physics of the Peierls Hubbard model is described by the spin-Peierls Heisenberg Hamiltonian with a staggered exchange interaction and a dimerization induced spin gap .

### B Correlation functions

The important question remains about the nature of the insulating phase of the IHM for $`U>U_c`$. To further analyze the BI and CI phases below and above $`U_c`$, we have evaluated site- and bond-charge distribution functions as well as spin-spin correlation functions.

In Fig. 4 we show the charge distribution $`0|(n(r)1)|0`$ ($`n(r)=n_{r,}+n_{r,}`$) in the ground state $`|0`$ for the IHM at $`U=0.8t<U_c`$ and $`U=4t>U_c`$ for a $`L=32`$ chain. The alternating pattern in the density distribution is well pronounced not only in the BI but also in the CI phase far beyond the critical point at $`UU_c`$. For the $`L=32`$ chain the CDW is well established at distances $`lL/2`$ even at $`U=4t`$ and its amplitude remains almost unchanged in the finite-size scaling analysis (see upper inset in Fig. 5). The main plot in Fig. 5 shows that the staggered component of the charge density decreases smoothly with increasing $`U`$. Close above the transition point near $`U=2t`$ one observes an anomaly, i.e. a slight decrease of the CDW amplitude accompanied by a change in curvature; this anomaly can be clearly identified in the enlargement shown in the lower inset. The anomaly will find a natural explanation in the discussion below. Our numerical data show that the alternating pattern in the electron density distribution in the IHM remains for arbitrary finite $`U`$. Thus, the ionic potential enforces long range CDW order for all interaction strengths.

Fig. 6 shows the DMRG results for the spin-spin correlation function $`0|S^z(L/2)S^z(L/2+r)|0`$ for the IHM at $`U=0.8t<U_c`$ and $`U=4t>U_c`$ in comparison with the Hubbard chain. In the BI phase at $`U=0.8t`$ the SDW correlations are quickly suppressed after a few lattice spacings. At $`U=4t`$ the amplitude of the SDW correlations in the CI phase of the IHM is slightly reduced in comparison to the Hubbard model at the same value of $`U`$. However, the large distance behaviors of the spin correlations in the CI phase of the IHM and the MI phase of the Hubbard model are quite similar and become almost indistinguishable (see the data in the inset for $`U=4t`$). On the other hand, the long distance behavior of the spin-spin correlation function at $`U=2.1t`$ and $`\mathrm{\Delta }=0.5t`$, i.e. close above the transition, manifestly supports the finiteness of the spin gap (see inset). This may be viewed as an indication for existance of two different phases above $`U_c`$. But we cautiously point out that it is hard to judge on the persistence or vanishing of $`\mathrm{\Delta }_S`$ far above $`U_c`$ in the CI phase from the finite chain spin correlators alone.

To address the BOW ordering tendencies in the CI phase we have calculated the ground-state distribution of the bond-charge density Eq. (23). Fig. 7 shows the results of the DMRG calculations for the $`L=32`$ IHM and the Hubbard chain at $`U=0.8t`$ and $`U=2.6t`$. The boundary effect for an open chain is strong and leads to a modulation of the bond density already for the pure Hubbard model. Interestingly, the same behavior was also observed previously for the bond expectation value $`𝐒_i𝐒_{i+1}`$ in open antiferromagnetic Heisenberg spin-$`1/2`$ chains . The detailed comparison with the Hubbard chain in Fig. 7 reveals, that the ionic potential leads to a reduction of the bond-density oscillations at $`U<U_c`$, while in the CI phase at $`U>U_c`$ their amplitude slightly increases. The enhancement of BOW correlations above $`U_c`$ must simultaneously weaken the CDW amplitude. This naturally explains the slight downward curvature near $`U_c`$ in the staggered charge density component shown in Fig. 5.

In order to explore the possibility towards true long range BOW ordering in the IHM above $`U_c`$ we plot in Fig. 8 the staggered bond-density versus $`U`$ in the center of long, open chains with $`L=256`$ and $`L=512`$. In the BI phase this quantity is essentially zero. At the transition point the staggered component of the bond density increases rapidly, and on further increasing $`U`$ it starts to decrease smoothly. However, the staggered bond density remains finite for any $`U>U_c`$ on these long but finite chains and vanishes only in the limit $`U\mathrm{}`$.

Naturally it is necessary to perform a finite size scaling analysis for this quantity. Fig. 9 shows its chain length scaling behavior for $`U=0`$, $`U=2.1t`$, and $`U=4t`$. Two conclusions can be drawn from these results: In the absence of the interaction the staggered bond density clearly extrapolates to zero - as expected for a conventional BI. In the CI phase close above the transition the upward curvature of the staggered bond density vs. $`1/L`$ points to a finite value in the infinite chain length limit, i.e. long range BOW order. However, the scaling behavior changes in the CI phase far above the transition point. Interestingly, the scaling behavior in the IHM far above $`U_c`$ starts to reseamble the results for the pure Hubbard model, for which the staggered bond density has to vanish in the thermodynamic limit. The qualitative change in the scaling behavior of the staggered bond density may again be viewed indicative for a possible second phase transition. Summarizing the results for the SDW and BOW correlations we conclude that there is evidence for two phases for $`U>U_c`$. Close above $`U_c`$ long range BOW order develops in the ground state of IHM while far above the transition point the coorrelation functions become almost identical to those of the Hubbard model. Yet, long range CDW order exists for all $`U`$. A precise location of a second transition point is however not possible from the currently awailable DMRG data.

## V Conclusions

From the finite chain DMRG studies and finite size scaling analyses we draw the following conclusions for the ground-state phase diagram of the IHM: the ionic potential leads to long range CDW order for all interaction strengths. The data resolve one transition point from the BI to the CI phase. Remarkably, at the transition $`\mathrm{\Delta }_C=\mathrm{\Delta }_S`$ and both remain finite. Close above the transition, i.e. for $`0<(U/U_c)11`$, we identify a clear signal for long range BOW order. With increasing $`U`$ above $`U_c`$ the finite size scaling behavior of the staggered bond density and spin-spin correlation function changes qualitatively and approaches the scaling behavior of the Hubbard model. With the current chain length and DMRG accuracy limitations it is not possible to precisely identify and locate the second transition point. The phase with BOW order necessarily has a finite spin excitation gap. If BOW seizes to exist above a second critical value of U, the spin gap has to vanish simultaneously leading to identical large distance decays of SDW and BOW correlations.

The insulator-insulator transition at $`U_c(\mathrm{\Delta })`$ on finite periodic chains results from a ground-state level crossing of the two site-parity sectors. The optical excitation gap therefore has to vanish at $`U_c`$; remarkably, the DMRG data reveal that at the critical point $`\mathrm{\Delta }_C`$ ans $`\mathrm{\Delta }_S`$ remain both finite and equal. This itself clearly indicates the existanse of a CI phase with $`\mathrm{\Delta }_C>\mathrm{\Delta }_S>0`$ which originates from the appearance of long range staggered bond density order. The distinction between the optical and the charge gap is therefore of key importance for the structure of the insulating phases and the phase transitions of the IHM. The investigation of the optical conductivity in the critical region is therefore a demanding task for future work on the complex physics of the ionic Hubbard model.

###### Acknowledgements.

We thank R. Noack, B. Normand, H. Fehske, and T. Vekua for helpful discussions. This work was supported by the Deutsche Forschungsgemeinschaft (DFG) through SP 1073. APK also acknowledges support through Sonderforschungsbereich 484 of the DFG. GIJ and MS acknowledge support by the SCOPES grant N 7GEPJ62379.


format: text

text:
# Stability and electronic structure of the (1×1) SrTiO3(110) polar surfaces by first principles calculations

## I Introduction

Strontium titanate is a material that has many potential applications. Its high dielectric constant makes it a candidate to replace silicon dioxide in some nano-electronic devices. SrTiO3 has been used in photo-electrolysis and as a substrate for the growth of high $`T_c`$ superconductors as well as of other thin oxide films Noguera (1996). Like other perovskite compounds, SrTiO3 shows a rich physical behavior. It undergoes an anti-ferro-distortive (AFD) transition at $`T<105K`$: the cubic structure turns into a tetragonal phase where the neighboring TiO6 octahedra are tilted with opposite angles around a  direction, doubling the unit cell. Moreover, SrTiO3 is at the boundary of a ferroelectric (FE) transition Itoh et al. (1999); Ravikumar et al. (1995); Pertsev et al. (2000) but remains paraelectric for all temperatures Zhong and Vanderbilt (1996).

While the (100) surfaces of SrTiO3 have been extensively studied, both theoretically Cheng et al. (2000); Padilla and Vanderbilt (1998); Heifets et al. (2001), and experimentally Jiang and Zegenhagen (1999); Kubo and Nozoye (2001); Castell (2002), the (110) terminations are much less known. Such a scarcity is likely due to the polar character of the (110) orientation. The sequence of atomic layers of O2 and SrTiO stoichiometry implies a monotonic raise of the microscopic electric field, which has to be compensated either through a modification of the surface composition - which leads to non-stoichiometric terminations - or by an anomalous filling of the surface electronic states - which must imply crucial variations of the electronic structure of the surfaces that should be in principle detectable by experiments.

Nevertheless, {110} terminations of strontium titanate has been observed quite often, but the effort towards the precise characterization of their atomic-scale morphology and the detailed study of the corresponding electronic structures started only in the last decade. For instance, it is worth noting the investigations by atomic force microscopy Szot and Speier (1999), scanning tunnelling microscopy (STM), Auger spectroscopy and low energy electron diffraction(LEED) measurements Brunen and Zegenhagen (1997), as well as ultraviolet photoemission and x-ray photoemission spectroscopies (UPS and XPS, respectively), coupled to LEED Aiura et al. (1994). In the latter works the authors made a considerable effort to characterize the modifications of the electronic structure that are connected to different preparation conditions of the surface through STM and scanning tunnelling spectroscopy (STS) Bando et al. (1995), angle resolved photoemission Haruyama et al. (1996) and conduction measurements Bando et al. (2001). Not all these studies provide a unifying picture of the surface, which indeed shows a great sensitivity to the thermodynamic conditions, mainly temperature and oxygen partial pressure.

Motivated by that experimental work, some theoretical studies of SrTiO3(110) recently appeared. The optimized surface geometry was obtained through interatomic forces that were derived either by an empirical shell-model Heifets et al. (2000), or by more accurate treatments that are based on the Hartree-Fock approximation, still describing the hamiltonian matrix on semi-empirical grounds Pojani et al. (1999); Stashans and Serrano (2002). The approximate nature of those approaches allowed some trends to be deduced, but it did not permit in most cases to draw quantitative conclusions that could be compared to experiments. On the other hand, first-principle approaches that are based on the Density Functional Theory (DFT) have been widely applied to the study of oxide surfaces in recent years. More specifically to perovskite surfaces, such an approach has been showed to be reliable and complementary to experiments Padilla and Vanderbilt (1998). It is the aim of the present paper to provide an application of state-of-the-art first-principle methods to the (110) polar termination of SrTiO3, which is the first one to our knowledge.

In particular, we focus on two issues: (i) Which is the mechanism at work to cancel the macroscopic polarization? and (ii) Is the thermodynamic stability of the (110) termination comparable to that of the (100) cleavage face? Regarding (i), while in the wide-gap rock-salt compounds MgO and NaCl the non-stoichiometric reconstructions such as the octopolar one are clearly favored over polarity compensation by an anomalous filling of surface states, it has been recently suggested for the less ionic ZnO$`(0001)`$ and ZnO$`(000\overline{1})`$ surfaces Wander et al. (2001); Carlsson (2001) that the latter mechanism may be at work. The coexistence of the rather ionic Sr-O bonds with the more covalent Ti-O ones, as well as the fairly large dielectric constant of SrTiO3 might situate strontium titanate at the borderline between the two previous cases. Regarding (ii), the quite numerous observations of {110} $`(1\times 1)`$ terminations suggest that these surfaces may be obtained in appropriate thermodynamic conditions. However, the ternary nature of strontium titanate, as well as the lack of precise experimental indications on the surface composition, make the first-principle modelling a rather delicate task. Therefore, in this work, we restrict ourselves to several terminations, stoichiometric or not, of $`(1\times 1)`$ reconstructions, completing a previous preliminary report Bottin et al. (2003): the two stoichiometric O2 and SrTiO faces, and the three non-stoichiometric O, TiO and Sr terminations. Their atomic and electronic structures are computed and discussed in Section III. Moreover, their relative stabilities are compared as a function of the chemical environment (see Sections IV and V).

## II Computational Method

The calculations are carried out within the density functional theory (DFT) Hohenberg and Kohn (1964). The exchange and correlation energy is treated via the local density approximation (LDA) using the Perdew-Wang parametrization Perdew and Wang (1992). The Kohn-Sham orbitals are expanded on a plane wave (PW) basis set and an energy cutoff of 30 Hartree (Ha) is employed. For all calculations, we use the ABINIT computer code Gonze et al. (2002).

In conjunction with the use of a PW basis set, we adopt soft norm-conserving pseudopotentials that are generated following the Troullier-Martins scheme Troullier and Martins (1991) in the Kleinman-Bylander form Kleinman and Bylander (1982) in order to avoid to take into account the inner atomic electrons in the self-consistent cycle explicitly. Along the pseudopotential generation process, atom reference states 3$`d^1`$4$`s^{0.5}`$4$`p^{0.25}`$ for Ti, 4$`s^2`$4$`p^6`$5$`p^1`$ for Sr and 2$`s^2`$2$`p^4`$ for O are used. The pseudization radii are $`R_c(s)`$=2.45 bohr, $`R_c(p)`$=2.55 bohr and $`R_c(d)`$=2.25 bohr for Ti, $`R_c(s)`$2.0 bohr, $`R_c(p)`$=2.5 bohr and $`R_c(d)`$=2.2 bohr for Sr and $`R_c(s)=R_c(d)`$=1.4 bohr and $`R_c(p)`$=1.75 bohr for the O, with $`s`$, $`p`$ and $`d`$ being the different channels. The Sr 4$`s`$ and 4$`p`$ semi-core electrons are treated as valence electrons in the self-consistent procedure, (which we refer to in the following as the small-core approximation) whereas the Ti 3$`s`$ and 3$`p`$ semi-core states are frozen (large-core approximation). Since the latter states have an appreciable superposition with the $`4s`$ and $`3d`$ radial orbitals of Ti, in the exchange and correlation functional we take into account the non-linear core correction Fuchs and Scheffler (1999) between the radial density corresponding to the frozen $`3s`$ and $`3p`$ Ti orbitals and the valence electron density.

The previously described computational scheme is tested on the cubic phase of bulk SrTiO3. The lattice parameter $`a_0`$, the bulk modulus $`B_0`$ and the cohesive energy $`E_{coh}`$ are very close to the experimental data (see Table 1). The large-core approximation for Ti is also checked against the small-core one, and the former is adopted in all slab calculations. The gap between the valence and conduction band computed through the difference of the Kohn-Sham eigen-energies (2.1 eV) is smaller than the experimental value (3.2 eV) Cardona (1965), which is the common case for the DFT Godby et al. (1988).

The calculated formation energies of bulk Ti and Sr from the isolated atomic species are: $`E_{\mathrm{Ti}}^f=4.58eV`$ and $`E_{\mathrm{Sr}}^f=1.73eV`$, which differ by -6$`\%`$ and +2$`\%`$ from the experimental values, respectively. On the other hand, the O2 formation energy ($`E_{\mathrm{O}_2}^f=7.58eV`$) is about 45$`\%`$ greater than the experimental value, like in previous LDA calculations Jones and Gunnarsson (1989).

The surfaces are described in the framework of the slab model. The calculations are carried out by sampling the irreducible Brillouin zone by a (4,4,2) Monkhorst-Pack mesh Monkhorst and Park (1976). By using a (5,5,2) mesh, the computed total energies decrease by less than 10 meV for the insulating terminations and 20 meV for the metallic ones. In those cases, when an effective Fermi-surface smearing is employed to get a better convergence Mermin (1965), a careful procedure is used to extrapolate convergent energies in the limit $`T0`$. The stable surface configurations are obtained by minimizing the Hellmann-Feymann forces in the Born-Oppenheimer approximation, generally starting from the geometry of the ideal unrelaxed surfaces.

In order to avoid spurious interaction between periodic slabs via dipole-dipole interactions Neugebauer and Scheffler (1992); Bengtsson (1999), the compositions of the adopted slabs are symmetric upon inversion along the surface normal. Well converged slab total energies can be routinely obtained for typical vacuum widths ranging between 10 and 12 Å. On the other hand, the determination of the minimal number of atomic layers that is allowed needs some care. For instance, a nine-layers slab is sufficient for the two non-stoichiometric (110)-Sr and (110)-TiO terminations, whereas 11 layers are at least needed for the two (110)-O2 and (110)-SrTiO stoichiometric terminations. In the case of the O-terminated surfaces, a very careful choice of the bulk reference energy is also needed, as detailed in the following. To this purpose, we made calculations up to 15-layers thick slabs. As a general rule, the slab is considered to be thick enough when the difference between the total energies of two slabs differing by a SrTiO3 unit and the bulk total-energy is smaller than 1 mHa 11endnote: 1 The uncertainty associated to this value, because of the number of bulk layers considered and the convergence with respect to the number of k points, is smaller than the thickness convergence one.. This condition cancels out the problem of nonconvergent surface energies Boettger (1994); Fiorentini and Methfessel (1996) and makes it possible to obtain a fair bulk reference energy for each termination.

Even if our calculations are carried out at 0 K, we impose symmetry constraints to our slabs to simulate the cubic phase that is stable at room temperature. A special case is represented by the non-stoichiometric (110)-O termination, since the mirror symmetry along the $`[\overline{1}10]`$ direction is lost. As it will be detailed in the section devoted to such a termination, the bulk energy should be computed with reference to an anti-ferro-distortive (AFD) phase that is simulated in a tetragonal supercell analogous to that used for the slab. Other authors gave details about the change of the AFD distortion with respect to the bulk structure Sai and Vanderbilt (2000). Our computed total-energy differences between the cubic and the AFD phases is about 2 mHa. Neglecting such a tiny quantity may appreciably bias the numerical extrapolation of the surface energy when a linear scheme is used. To summarize, the uncertainty on the surface energies $`E_{cl}`$ within our slab approach does not exceed $`10^2`$ J/m2, which is usually much smaller than the computed surface energy differences for the various terminations.

We also focus on the total electron distribution in the various surface configurations. An estimate of the electron sharing between the O anions and the Sr and Ti cations is obtained by means of the Bader’s topological analysis of the electron density, which corresponds to a partition of the total charge in atomic bassins. Such a procedure is independent of the basis set used Bader (1991), and has been recently used Noguera et al. (2002) and improved within the ABINIT package. The precision of the atomic charge integration is about $`5\times 10^3`$ electrons.

## III Electronic and atomic structure

In the (110) orientation of SrTiO3, a stacking sequence of atomic layers …O2-SrTiO-O2-SrTiO… is provided. If we consider that the ionic charge of O, Ti and Sr are $`Q_\mathrm{O}=2`$, $`Q_{\mathrm{Ti}}=+4`$ and $`Q_{\mathrm{Sr}}=+2`$ respectively, then the O2 plane bears a formal charge $`Q_{\mathrm{O}_2}=4`$ and the SrTiO plane $`Q_{\mathrm{SrTiO}}=+4`$ per 2D unit cell. According to the criterion for polarity compensation Noguera (2000); Noguera et al. (2002), the formal surface charges of the various SrTiO3(110) faces have to be equal to half the bulk value, i.e., $`Q_{surf}=\pm 2`$. Therefore, we consider two main classes of $`(1\times 1)`$ terminations: on one side, the so called stoichiometric terminations, since their compositions reflects the bulk stacking: the (110)-SrTiO and -O2 ones, for which an anomalous filling of surface states is expected. On the other side, we study the (110)-TiO, -Sr and -O terminations, for which the stoichiometry changes can in principle provide the polarity compensation. For each termination, we describe in the following the detailed atomic and electronic structure.

### III.1 Stoichiometric terminations

#### III.1.1 The SrTiO termination.

In the stoichiometric (110)-SrTiO termination (see Figure 1), 6 Sr-O bonds and 2 Ti-O bonds are cut. Since the surface plane bears a formal charge $`Q_{\mathrm{SrTiO}}=+4`$ there must be an anomalous filling of surface states, in order to cancel out the macroscopic component of the slab dipole Noguera (2000). Figure 2 reports the computed band structure. One can note that the Fermi energy is located above the bottom of the conduction band. Some conduction-like states are thus filled, which bring the two additional electrons that are needed according to the electron counting rule for polarity compensation.

These states are localized in the outermost layers and are genuine surface states. The density of the filled conduction band state at the M point of the Brillouin zone that is displayed in figure 3, is roughly delocalized over three surface layers, which shows that the self-consistent charge redistribution is more complex than the model that is based upon the values of the formal ionic charges. A more thorough insight can be obtained by means of the Bader’s topological charge analysis, which is reported in Table 2. The layer charge q$`{}_{}{}^{I,III}{}_{\mathrm{SrTiO}}{}^{}`$ and q$`{}_{}{}^{II}{}_{\mathrm{O}_2}{}^{}`$ of the first and third SrTiO, and the second O2 outermost planes, respectively, are indeed strongly modified by the filling of the additional surface state and thus show marked changes from the computed charge for the inner bulk-like layers (q$`{}_{}{}^{Bulk}=\pm 2.50`$), which is roughly achieved only in the fourth layer. The polarity compensation criterion is fulfilled, since the sum of the charge of the four first layers is equal to -1.30, which is almost the half of the bulk layer charge qBulk.

Finally, we can see in table 2 that the difference between the surface layer charge q$`{}_{}{}^{I}{}_{\mathrm{SrTiO}}{}^{}`$ and the bulk layer charge qBulk is mainly localized on the Titanium, which is confirmed by the direct visualization of the filled state at the conduction band bottom (in the right part of figure 4). However, an effective decrease of the electronic kinetic energy may be obtained through the delocalization of this state over neighboring sites. Such a screening phenomenon, reminiscent to that found in Na-covered TiO2(110) surfaces Albaret et al. (2002), thus affects the Sr and O charges, too. The charge of the outermost Ti is thus reduced, which can be probed by XPS experiments through a surface Ti signal that should differ from the bulk Ti4+ formal oxidation state to a great extent.

As far as the atomic structure is concerned, we can see in table 3 that large relaxations happen in the four outermost layers. For instance, the large rumpling between Oxygen and Strontium that is equal to 0.73 Å in the surface layer, is still non negligible in the third layer (0.32 Å). The inter-planar distance between the outermost SrTiO layer and the O2 layer beneath contracts by 11 %, while that between the latter and the third SrTiO planes grows by +19%, showing a typical damped oscillation that makes the surface unlike layers closer two-by-two, which is accounted for by the general theory Desjonquères and Spanjaard (1993); Allan and Lannoo (1988). However, the surface Ti-O bondlengths $`d_{\mathrm{Ti}\mathrm{O}}`$ are only weakly modified: for Ti in the topmost layer, $`d_{\mathrm{Ti}\mathrm{O}}`$ is reduced of -3%, while it is expanded by 5% in the third layer. Such a behavior is due to the combined relaxation and rumpling that mainly rotate the Ti-O bond. At odds, the more ionic and weaker Sr-O bonds are much more affected: $`d_{\mathrm{Sr}^\mathrm{I}\mathrm{O}^{\mathrm{III}}}`$ is reduced by 13%, while $`d_{\mathrm{O}^\mathrm{I}\mathrm{Sr}^{\mathrm{III}}}`$ is expanded by 24%

#### III.1.2 The O2 termination.

The complementary stoichiometric (110)-O2 termination is also expected to undergo deep modifications of the electronic structure because of the polarity compensation, which requires two electrons to be removed from the surface. As displayed in the left part of Figure 5, the two surface oxygens move close to each other forming a bond 1.48 Å long. Such a bond-length is usually a signature of a peroxo group Cotton and Wilkinson (1980), which is formally denoted as $`\mathrm{O}_2^2`$ in the ionic limit. It would imply that the surface charge per unit cell is $`Q_{surf}=2`$ instead of -4 like in the bulk. Therefore, the formation of a peroxo group formally permits to fulfill the polarity compensation criterion.

Less formally, one can see that the formation of a peroxo bond is made possible by emptying the anti-bonding $`\sigma ^{}`$ molecular orbital (see the right panel of Figure 5). This also provides an effective mechanism to open a small gap at the surface (see figure 6). The Fermi level rests above the two anti-bonding states $`\pi ^{}`$ and does not cross the band structure . The bulk gap value of SrTiO3 (around 2 eV) is found between the valence band and the conduction band (i.e. excluding the O2 $`\pi ^{}`$ states). The two in-gap $`\pi ^{}`$ states are split since the two $`[\overline{1}10]`$ and  crystallographic directions are not equivalent, and are localized at the surface layer as displayed in figure 7.

Table 4 yields the results of a Bader’s topological analysis for this termination. The topological charges of the surface O ($`\mathrm{O}^I`$) are sensitively reduced, consistently with the formation of the peroxo bond. We obtain a surface charge almost equal to half the bulk one. In this case, the charge modification is essentially restricted to the outermost O atoms, as a consequence of the covalent and localized nature of the peroxo bond.

As far as the atomic structure is concerned, important modifications are present on the surface. The Oxygens tilt at the surface and reduce their angle with respect to the  direction by strongly increasing the bond-length with the Titanium on second layer ($`d_{\mathrm{O}^\mathrm{I}\mathrm{Ti}^{\mathrm{II}}}`$) by 19 %. As a consequence,there is a big inward relaxation (-20$`\%`$ – see also table 5) of the peroxo bond with respect to the second layer. On the other hand, $`d_{\mathrm{O}^\mathrm{I}\mathrm{Sr}^{\mathrm{II}}}`$ is remarkably reduced (-13 %).

We point out that a few surface configurations were obtained for the (110)-O2 termination through the minimization procedure. The one reported in Fig.5 has the lowest total energy. Another configuration with the peroxo group on top of the Titanium is found, which is very slightly higher in energy (5.$`10^3J/m^2`$) but still within the intrinsic precision of our calculations. These two configurations have very similar band structures and atomic charges, and can be considered as almost degenerated. It is worth noting that when starting the geometry optimization from the ideal, unrelaxed (110)-O2 termination, a third local minimum without peroxo groups at the surface is found. With respect to the peroxo configurations, it has an open-shell structure and a higher surface energy of about 0.7 J/m2.

### III.2 Non-stoichiometric terminations

The (110)-TiO, -O and -Sr terminations may be ideally obtained from the stoichiometric surface by adsorption or desorption of Sr and O atoms. They bear a formal charge $`Q_{\mathrm{TiO}}=+2`$, $`Q_\mathrm{O}=2`$ and $`Q_{\mathrm{Sr}}=+2`$ in the ionic limit, respectively, which corresponds to half the bulk layer charge $`Q_{bulk}=\pm 4`$. Therefore, they are all compensated and no anomalous filling of surface states is in principle necessary. In order to demonstrate this assumption, the computed density of states of those three non-stoichiometric terminations are drawn in Fig.8. They are all insulating, with electronic structures qualitatively similar to the bulk.

#### III.2.1 The Sr and TiO terminations.

The (110)-Sr termination is obtained from the O2 one by adsorbing a Sr atom per unit cell whereas the (110)-TiO one is obtained from the (110)-SrTiO termination by removing a row of surface Sr along the  direction. These two terminations are displayed in figure 9.

These two (110) terminations are representative of more open surfaces. The coordination numbers of the surface atoms are reduced with respect to the stoichiometric (110) $`(1\times 1)`$ terminations as well as to the non-polar (100) surfaces. Let us imagine to cleave a SrTiO3 slab along different orientations: two complementary terminations are obtained, for which a certain number of bonds are missing. Such cut bonds are: 4 Sr-O and 1 Ti-O for the non-polar (100) orientation; 2 Ti-O and 6 Sr-O bonds for the (110)-SrTiO and (110)-O2 terminations, as well as for the (110)-O termination, which is self-complementary; 2 Ti-O and 8 Sr-O bonds for the (110)-TiO and (110)-Sr terminations. Moreover, at variance with the TiO2 and SrO (100) terminations and SrTiO, O2 and O (110) ones, the (110)-Sr and the (110)-TiO terminations are the only ones having an under coordinated atom even in the third layer starting from the surface (a Sr-O bond is cut). Therefore, it is not surprising that the remaining atoms on the (110)-Sr and (110)-TiO undergo large relaxations (see table 6), like the inward relaxation of the surface Sr on the Sr termination (-43$`\%`$) and the strong rumpling of the Ti and O atoms on the TiO termination (0.48$`\AA `$). Despite this under-coordination, relaxations are more quickly damped than for the stoichiometric terminations. For instance, the inter-planar relaxation between the third and the fourth layer is reduced by a factor two with respect to the stoichiometric (110)-SrTiO and (110)-O2 terminations.

Since the (110)-Sr and the (110)-TiO terminations are intrinsically compensated by stoichiometry, we may guess that the atomic surface charges differ little from the bulk, and that such an effect is mainly be due to the reduction of the surface atom coordination numbers. In fact, the Bader’s topological analysis (see Table 7) confirms the polarity compensation of the two terminations: the charge modification affects essentially the three outermost layers and their sum on these layers is equal to +1.21, which is almost half of the bulk value. Even if these changes are less intense than on the stoichiometric terminations, they extend rather deeply into the slabs, which confirms the correlation with the presence of under-coordinated atoms and the interplay between the atomic relaxations and the electron redistribution.

#### III.2.2 The O termination.

Such a termination may be ideally obtained in two different ways: by removing an Oxygen atom from a O2 termination or by adding an Oxygen atom on a SrTiO termination. The presence of a single O atom on the surface breaks the mirror symmetry along the $`[\overline{1}10]`$ direction, at odds with all $`(1\times 1)`$ terminations that were previously considered. Consequently, among all the (110) terminations, the convergence of the atomic structure and the surface energy of this one is the most delicate. We can see in figure 10 that the broken surface symmetry permits a distortion and its propagation into the slab. It consists in an alternating rotation of the octahedra along a  direction, which is known as an anti-ferro-distortive (AFD) instability. We point out that such a abnormal deep propagation of the AFD instability from the surface into the bulk is related to the fact that the low-symmetry phase is stable at T=0 K, the temperature for which the calculations are performed. This propagation may be hindered by considering a 15 layer slab with the five inner layers frozen. Even in that case, the surface atomic structure is practically undistinguishable from that obtained on top of an AFD bulk.

The surface O relaxes inwards (-33$`\%`$) and laterally toward the two Strontium underneath (see table 8). In the inner part of the slab, no large inter-planar relaxation happens. However, a quite strong rumpling affects the inner O2 layers, which corresponds to the AFD distortion and is not damped as a function of the slab thickness. Moreover, there are large lateral atomic displacements within the SrTiO layers that are denoted with $`\delta _{\mathrm{TiO}_2}`$ (see table 8) and decrease quite rapidly going into the bulk. As a consequence of this distortion, the O-Ti-O alignement is broken along the  direction.
The topological charges of the surface oxygen’s are very similar to the bulk ones, from which they differ by 10% at most.

## IV Thermodynamic stability

Because of the different chemical nature of the various $`(1\times 1)`$ (110) terminations of SrTiO3, their stability must be discussed as a function of the actual surface composition Wang et al. (2000, 1998); Padilla and Vanderbilt (1998); Pojani et al. (1999). To this purpose, we make use of two distinct physical quantities: the first one corresponds to the energy that is necessary to split the crystal in two parts and create two complementary terminations, which we refer to as the cleavage energy $`E_{cl}`$. The second one is the surface grand-potential $`\mathrm{\Omega }_s^i`$, which is a measure of the excess energy of a semi-infinite crystal exposing a termination with a given composition $`i`$ in contact with matter reservoirs.

### IV.1 Cleavage energy

When a stoichiometric SrTiO3 slab is ideally cut and the two parts are put apart, two complementary surface terminations are created. In our case, they are the (110)-SrTiO and the (110)-O2 terminations on one hand, the (110)-TiO and the (110)-Sr ones on the other. The (110)-O termination can be considered as self-complementary. The respective cleavage energies $`E_{cl}^{(\mathrm{O}_2+\mathrm{SrTiO})}`$ and $`E_{cl}^{(\mathrm{TiO}+\mathrm{Sr})}`$ can be obtained from the total energies computed for the symmetric slabs through the following equations:

$`E_{cl}^{(\mathrm{TiO}+\mathrm{Sr})}`$ $`=`$ $`{\displaystyle \frac{1}{2S}}(E_{slab}^{\mathrm{TiO}}+E_{slab}^{\mathrm{Sr}}nE_{bulk})`$ (1)
$`E_{cl}^{(\mathrm{O}_2+\mathrm{SrTiO})}`$ $`=`$ $`{\displaystyle \frac{1}{2S}}(E_{slab}^{\mathrm{O2}}+E_{slab}^{\mathrm{SrTiO}}nE_{bulk})`$ (2)
where $`E_{slab}^i`$ the total-energy of the symmetric slab with the $`i`$ termination, $`n`$ the total number of bulk formula units in the two slabs, $`S`$ the surface area and Ebulk the bulk energy per formula unit in the cubic structure.

Since the (110)-O termination is self-complementary and shows an AFD distortion, its cleavage energy $`E_{cl}^{(O+O)}`$ is calculated by using a unit factor instead of one half and a different reference energy E$`{}_{}{}^{AFD}{}_{bulk}{}^{}`$, which reads:

$$E_{cl}^{\mathrm{O}+\mathrm{O}}=\frac{1}{S}(E_{slab}^\mathrm{O}nE_{bulk}^{AFD})$$
(3)
As previously pointed out in Section II, a different choice for the reference energy would result in a ill-defined cleavage energy. This is exemplified in a recent calculation, in which an apparently diverging cleavage energy was obtained for the (110)-O termination Heifets et al. (2000) by using the cubic, undistorted reference energy instead of $`E_{bulk}^{AFD}`$. Our numerical results are summarized in table 9. The O termination has the lowest cleavage energy (2.54 J/m2) whereas the O2 and SrTiO the highest one (6.52 J/m2). We also note that the three non-stoichiometric terminations have lower cleavage energies than the two stoichiometric ones. Consequently, the polarity compensation that is achieved through the modification of the surface stoichiometry seems to be more effective than that by the anomalous filling of the surface states, as far as the energetics is concerned.

### IV.2 The surface grand-potential.

In order to distinguish the contribution of each termination to the cleavage energy, we compute its surface grand potential, which implies a contact with matter reservoirs. Many authors have recently used this method with success, as well for binary Wang et al. (2000, 1998) as for ternary compounds Padilla and Vanderbilt (1998); Pojani et al. (1999). We introduce the chemical potential $`\mu _{\mathrm{Ti}}`$, $`\mu _{\mathrm{Sr}}`$ and $`\mu _\mathrm{O}`$ of the Ti, Sr and O atomic species, respectively. The surface Grand Potential per area unit $`\mathrm{\Omega }^i`$ of the $`i`$ termination reads 22endnote: 2The -TS term is neglected for typical temperatures T (S is the entropy of the system):

$$\mathrm{\Omega }^i=\frac{1}{2S}[E_{slab}^iN_{\mathrm{Ti}}\mu _{\mathrm{Ti}}N_{\mathrm{Sr}}\mu _{\mathrm{Sr}}N_\mathrm{O}\mu _\mathrm{O}]$$
(4)
with $`N_{\mathrm{Ti}}`$, $`N_{\mathrm{Sr}}`$ and $`N_\mathrm{O}`$ the number of Ti, Sr and O atoms in the slab and the half factor corresponds to the surface Grand potential by termination. The chemical potential $`\mu _{\mathrm{SrTiO}_3}`$ of a condensed and stoichiometric phase of strontium titanate is written as a sum of three terms representing the chemical potential of each species within the crystal:

$$\mu _{\mathrm{SrTiO}_3}=\mu _{\mathrm{Sr}}+\mu _{\mathrm{Ti}}+3\mu _\mathrm{O}$$
(5)
Since the surface is in equilibrium with the bulk SrTiO3, we have $`\mu _{\mathrm{SrTiO}_3}=E_{bulk}`$. If we replace equation 5 in equation 4, we can eliminate the $`\mu _{\mathrm{Ti}}`$ and $`\mu _{\mathrm{SrTiO}_3}`$ variables in the surface Grand Potential and obtain:

$$\mathrm{\Omega }_s^i=\frac{1}{2S}[E_{slab}^iN_{\mathrm{Ti}}E_{bulk}\mu _\mathrm{O}(N_\mathrm{O}3N_{\mathrm{Ti}})\mu _{\mathrm{Sr}}(N_{\mathrm{Sr}}N_{\mathrm{Ti}})]$$
(6)
Relying upon equation 6, one can deduce, for each termination, the range of the accessible values of $`\mathrm{\Omega }_s^i`$ if the minimum and maximum values of the O and Sr chemical potentials are known (see appendix).

If we introduce the variation of the chemical potentials with respect to those computed for the reference phases ($`\mathrm{\Delta }\mu _\mathrm{O}=\mu _\mathrm{O}\frac{E_{\mathrm{O}_2}^{mol}}{2}`$ and $`\mathrm{\Delta }\mu _{\mathrm{Sr}}=\mu _{\mathrm{Sr}}E_{\mathrm{Sr}}^{bulk}`$, respectively) in equation 6 we obtain:

$`\mathrm{\Omega }_s^i`$ $`=`$ $`\varphi _i{\displaystyle \frac{1}{2S}}[\mathrm{\Delta }\mu _\mathrm{O}(N_\mathrm{O}3N_{\mathrm{Ti}})\mathrm{\Delta }\mu _{\mathrm{Sr}}(N_{\mathrm{Sr}}N_{\mathrm{Ti}})]`$
$`\text{with }\varphi _i`$ $`=`$ $`{\displaystyle \frac{1}{2S}}[E_{slab}^iN_{\mathrm{Ti}}E_{bulk}`$ (7)
$`{\displaystyle \frac{E_{\mathrm{O}_2}^{mol}}{2}}(N_\mathrm{O}3N_{\mathrm{Ti}})E_{\mathrm{Sr}}^{bulk}(N_{\mathrm{Sr}}N_{\mathrm{Ti}})]`$
$`\varphi _i`$ measures the stability of the surface with respect to bulk SrTiO3, gaseous molecular oxygen and metallic Sr. In Table 10 our ab initio results for $`\varphi _i`$, are listed.

For two complementary terminations (TiO-Sr and SrTiO-O2) the sum of their surface Grand potential is independent of the chemical potential and corresponds to their cleavage energy (see Table 9).

The derivation of the upper and lower bounds of $`\mathrm{\Delta }\mu _\mathrm{O}`$ and $`\mathrm{\Delta }\mu _{\mathrm{Sr}}`$ is detailed in appendix A. Within the allowed region, we show in figure 11 (left panel) the $`(1\times 1)`$ (110) terminations having the lowest surface grand potentials, which provides the stability diagram of SrTiO3(110) $`(1\times 1)`$ surfaces in an O and Sr external environment.

First of all, according to Fig.11 (left panel), our calculations predict that only four out of the five possible terminations may be obtained. Indeed, the (110)-O2 termination cannot be stabilized, even in very O-rich chemical environments. The Sr termination is the most stable one in O- and Sr-rich environments, as its complementary TiO face is in O- and Sr-poor conditions. The O termination shows a stability domain in moderate O- and Sr-environment. Finally, the stoichiometric and open-shell (110)-SrTiO termination happens to be stable in a small domain corresponding to O-poor and Sr-rich conditions. In order to dicuss the existence of this small domain (at least within the theory) with respect to the precision of our calculations, we point out that (see the right panel of Fig. 11):

* At the O-poor and simultaneously Sr-rich zone boundary (i.e.: $`\mathrm{\Delta }\mu _O`$=-5.46 eV and $`\mathrm{\Delta }\mu _Sr`$= 0 eV) the difference between Sr and SrTiO surface Grand Potentials is equal to 0.13 J/m2, which is higher than the estimated precision of our calculations, which is of the order of 0.01 J/m2.
* Since the area of the ($`\mathrm{\Delta }\mu _O`$, $`\mathrm{\Delta }\mu _Sr`$) domain is slightly underestimated by our first principle calculations (see appendix), the consideration of the experimental boundaries would enlarge the stability domain of the SrTiO termination.

## V Discussion

Among the questions that we raise in the introduction, in the following we discuss three issues that are especially noteworthy.

### V.1 A stable polar stoichiometric termination: (110)-SrTiO.

We have shown in the previous section that the $`(1\times 1)`$ stoichiometric (110)-SrTiO termination can be stabilized in a O-poor and Sr-rich environment. Even if its stability domain comes out to be rather small, this is the second case (to our knowledge) together with the Zn and O terminations of ZnO(0001) Wander et al. (2001) of a stable stoichiometric polar oxide surface with an open-shell electronic structure. However, a recent STM investigation Dulub et al. (2003) questioned the proposed stoichiometric morphologies for ZnO(0001) and provides different, non stoichiometric, structural models.

The a polarity compensation through anomalous filling of surface states is expected on the basis of theoretical arguments, but it is actually not often encountered Noguera (2000). Indeed, the cleavage energies of stoichiometric polar terminations in rock-salt structures, such as the $`(1\times 1)`$ unreconstructed (111) faces of MgO, are generally much higher than those of reconstructed, non-stoichiometric polar terminations. As a consequence, the surface grand potentials of these polar stoichiometric terminations result higher than those of the non stoichiometric ones even when the chemical potential dependent terms represent a negative contribution. In this respect, the peculiar behavior of SrTiO3 may be due to the presence of Ti-O covalent bonds, in conjunction with a not too large fundamental gap. Indeed, the energy increase that is due to the anomalous filling of surface states needed for polarity compensation may be effectively lessen by atomic relaxation and electronic screening effects. The latter ones show up through a non negligible charge transfer affecting some surface layers as supported by the analysis of the topological Bader’s charge of the SrTiO termination (Table 2): all Ti and O charges belonging to the three topmost layers are modified. Such a screening mechanism is less costly than a drastic charge reduction on only one or two surface atoms, which is essentially the case for MgO(111) $`(1\times 1)`$ Finocchi .

### V.2 Comparison with the non polar SrTiO3(100) terminations.

The cleavage (100) orientation of SrTiO3 is non-polar and displays two different stoichiometric terminations: the (100)-TiO2 and the (100)-SrO. They are usually expected to be more stable than any polar face such as the (110) terminations. Indeed, within the same theoretical and computational framework, we calculated the cleavage energy of the $`(1\times 1)`$ (100)-TiO2 and -SrO surfaces, which is not very much lower than those of the non-stoichiometric (110) terminations. Therefore, it is worth comparing the thermodynamic stability of the simulated $`(1\times 1)`$ (110) terminations with the two $`(1\times 1)`$ (100)-TiO2 and (100)-SrO faces. In figure 12 the domains of stability of the two (100) terminations and the five (110) terminations previously reported are gathered. Even if in the most common conditions, corresponding to moderate Sr and O chemical potentials, the (100) faces are favored, two (110) distinct $`(1\times 1)`$ terminations are predicted to be stable – the (110)-TiO in Sr-poor environments and the (110)-Sr in O-rich and Sr-rich conditions – the first of which shows a quite wide domain of thermodynamic stability. From the theoretical point of view, the (110)-TiO and the (110)-Sr terminations can be thought as exposing strongly relaxed {100} nano-facets (see Fig.4). It is therefore not completely surprising that their surface energies can become comparable to those of flat {100} terraces, if the formation energy of corners and edges is not too high. This is consistent with the experimental evidence, which shows that {110} orientations can be quite easily obtained for SrTiO3.

### V.3 Comparison with experimental results.

Using the first-principles stability diagrams displayed above, we can now discuss and propose a tentative explanation for the experimental measurements. As anticipated in the introduction, the SrTiO3(110) surface was produced and characterized by several groups, showing a great sensitivity to the preparation conditions. In particular, if the SrTiO3(110) is slightly annealed in ultrahigh vacuum (UHV), at 800$`{}_{}{}^{\mathrm{o}}\mathrm{C}`$ Bando et al. (1995) or at 960$`{}_{}{}^{\mathrm{o}}\mathrm{C}`$ during 2 hours Brunen and Zegenhagen (1997), the surface exhibits a (1$`\times `$1) LEED pattern. In contrast, a $`(n\times m)`$ periodicity is observed when the SrTiO3(110) surface is annealed at temperatures $`T1000{}_{}{}^{\mathrm{o}}\mathrm{C}`$. Therefore, one may guess that $`(1\times 1)`$ SrTiO3 (110) surfaces represent local minima of the possibly complicated free-energy landscape in a quite wide temperature range. Preliminary calculations on larger cell reconstructions Bottin et al.  show that some $`(1\times 1)`$ phases may remain stable with respect to other reconstructed surfaces. In the following, we mainly discuss the experiments performed on these $`(1\times 1)`$ phases.

As far as the composition of the surface is concerned, the authors of Ref. Brunen and Zegenhagen, 1997 get evidence through Auger measurements of an increasing surface concentration of Sr at the high-temperature annealed surface. Therefore, they conclude for a lack of Strontium for the mildly annealed (1$`\times `$1) surface, and propose a {100} micro-faceted (110)-TiO2 termination as the corresponding atomic-scale structural model. On the other hand, the STM study of Ref. Bando et al. (1995) shows on the (1$`\times `$1) terminations the presence of rather flat regions where tunnelling spectra have metallic character. These findings are consistent with typical Ti3+ and Ti2+ features in the XPS spectra, as well as a metallic surface state with a Ti 3d character as seen in UPS. On the basis of all those observations, and at odds with Ref. Brunen and Zegenhagen, 1997, they propose a SrTiO termination for the unreconstructed SrTiO3(110) surface, which can account for its flatness and metallic character dominated by Ti-like occupied states.

Relying on our calculations, we argue that these two apparently contradictory models are not necessarily incompatible. On one hand, the UHV Auger measurements may have been carried out in Sr-poor environments, whereas the UHV experimental conditions that have been used in Ref. Bando et al., 1995 imply an O-poor environment. Keeping in mind the great sensitivity of the SrTiO3(110) surface to the actual thermodynamic conditions, it is not surprising that its atomic-scale structure may show big variations. According to our calculations (see Fig.11), the (110)-TiO termination, which can be interpreted as a {100}-TiO2 micro-faceted surface in agreement with the model proposed by Brunen and Zegenhagen Brunen and Zegenhagen (1997), can be obtained in a wide domain corresponding to slightly Sr-poor and O-poor environments. If we assume that the temperature annealing is done in O-poor and, at the same time, rather Sr-rich conditions (whether one starts from the very beginning with a stoichiometric SrTiO termination without desorbing Sr atoms or the latter ones migrate from the bulk to surface domains of TiO composition remains an open question) a local (110)-SrTiO termination may be obtained as proposed by Bando and coworkers Bando et al. (1995). Indeed, the metallic character of the (110)-SrTiO termination that is due to a surface state mainly of Ti character agrees with the measured UPS spectra. Moreover, as we have previously pointed out for this termination, the anomalous filling of a Ti surface state that is needed for polarity compensation, should be associated to special Tin+ features in XPS, with $`n<4`$.

## VI Conclusion

The SrTiO3(110)(1$`\times `$1) surface was studied in the framework of ab initio calculations for the first time. The number of surface terminations in this ternary compound is larger than in binary compounds, which permits a remarkable variety for mechanisms of polarity compensation. An anomalous filling of the surfaces states takes place at the SrTiO termination with an open-shell electronic structure and at the O2 termination with in-gap states. If non-stoichiometric O, TiO and Sr terminations are considered, we obtain for each of them an insulating bulk-like electronic structure.

By calculating the surface grand potentials, we obtain four distinct $`(1\times 1)`$ stable (110) terminations as a function of the Sr and O chemical potentials. A quite large domain is found for the three non-stoichiometric terminations, especially for the TiO one. This competition between different terminations is also confirmed by the presence of a small domain of stability for the SrTiO termination in simultaneously Sr-rich and O-poor environments.

In conclusion, due to the scarcity of experimental investigations and the observed complexity of the physical and chemical behaviors of the SrTiO3(110) surface, the detailed determination of its atomic and electronic structure still remains an open question. However, our first-principles calculations suggest that even for the unreconstructed $`(1\times 1)`$ surfaces many distinct terminations are likely to appear, according to the precise experimental conditions. In particular, the great dispersion of the experimental results should be connected with such a sensitivity. In the light of our simulations, the available structural models that have been proposed to interpret the measurements seem to be reasonable within a restricted thermodynamic domain. We point out that coupling spectroscopic measurements and structure-sensitive techniques such as near-field microscopies and grazing x-ray diffraction, in very carefully controlled chemical environments may be crucial to get a comprehensive insight into SrTiO3 polar surfaces.

###### Acknowledgements.

We thank P. Casek for helping us in performing the Bader’s topological charge analysis. We also ackowledge financial support from the GDR CNRS ”Interface et Surface Sensible à la Structure” and from the Laboratoire de Physique des Solides in Orsay, France. The calculations were done on the IBM RS/6000 SP Power3 computer at IDRIS, under project No. 24089.

## Appendix A Boundary limits for the chemical potentials

In this appendix we derive the range of accessible values for $`\mu _{\mathrm{Sr}}`$ and $`\mu _O`$ in the stability diagram of the SrTiO3(110) orientation. The Oxygen, Titanium and Strontium atoms are assumed to form no condensate on the surface. Consequently, the chemical potential of each species must be lower than the energy of an atom in the stable phase of the considered species:

$`\mathrm{\Delta }\mu _O`$ $`=`$ $`\mu _O{\displaystyle \frac{E_{\mathrm{O}_2}^{mol}}{2}}<0`$ (8)
$`\mathrm{\Delta }\mu _{\mathrm{Sr}}`$ $`=`$ $`\mu _{\mathrm{Sr}}E_{\mathrm{Sr}}^{bulk}<0`$ (9)
$`\mathrm{\Delta }\mu _{\mathrm{Ti}}`$ $`=`$ $`\mu _{\mathrm{Ti}}E_{\mathrm{Ti}}^{bulk}<0`$ (10)
We have introduced in these inequations the relative values $`\mathrm{\Delta }\mu _O`$, $`\mathrm{\Delta }\mu _{\mathrm{Sr}}`$ and $`\mathrm{\Delta }\mu _{\mathrm{Ti}}`$ of the different chemical potentials with respect to $`E_{\mathrm{Sr}}^{bulk}`$, $`E_{\mathrm{Ti}}^{bulk}`$ and $`\frac{E_{\mathrm{O}_2}^{mol}}{2}`$, which are the energies of a Ti atom in the hcp bulk metal, of the Sr atom in the bulk cubic structure and of the O atom in the O2 molecule in the gas phase, respectively. The two first inequations 8 and 9 define the upper boundaries of the Oxygen and Strontium chemical potentials. By combining Eq. 10 and Eq. 5 in the main text, we obtain the following lower boundaries:

$`\mathrm{\Delta }\mu _{\mathrm{Sr}}+3\mathrm{\Delta }\mu _O`$ $`>`$ $`E_{\mathrm{SrTiO}_3}^f\text{ with}`$ (11)
$`E_{\mathrm{SrTiO}_3}^f`$ $`=`$ $`E_{\mathrm{SrTiO}_3}^{bulk}E_{\mathrm{Ti}}^{bulk}E_{\mathrm{Sr}}^{bulk}{\displaystyle \frac{3}{2}}E_{\mathrm{O}_2}^{mol}`$
$`E_{\mathrm{SrTiO}_3}^f`$ is the formation energy of SrTiO3 with respect to the Ti and Sr atoms in their bulk phases, and the O atom in the gas phase wich is positive defined. In order to compute easily this quantity and compare it with the experimental values, we rewrite $`E_{\mathrm{SrTiO}_3}^f`$ in the following way:

$$E_{\mathrm{SrTiO}_3}^f=E_{coh}E_{\mathrm{Sr}}^fE_{\mathrm{Ti}}^f\frac{3}{2}E_{\mathrm{O}_2}^f$$
(12)
with all these quantities defined in Sec. II. Our computed value is $`E_{\mathrm{SrTiO}_3}^f=16.38`$ eV to be compared to the value $`E_{\mathrm{SrTiO}_3}^f=17.40`$ eV that can be deduced by the experimental data.
In conclusion we show the stability of the Strontium Titanate diagram in a ($`\mathrm{\Delta }\mu _{\mathrm{Sr}}`$;$`\mathrm{\Delta }\mu _\mathrm{O}`$) plane in Fig. 13. As one can see, the ab initio calculations only slightly underestimate the size of the stability region.