Electronic correlations and crystal structure distortions in BaBiO₃

BaBiO₃ is a parent compound for high-T_c superconductors Ba_1−xK_xBiO₃ and Ba_1−xPb_xBiO₃. Properties of these compounds are the most intriguing when chemical composition is near the transition to the SC state. It would be promising to describe softening of phonon modes with doping as a result of ab initio calculation in order to support the phonon mechanism for superconductivity in these materials. The proper description of the electronic and crystal structure of pure BaBiO₃ has serious difficulties within density functional theory (DFT). A possible solution to this problem is to take into account strong electronic correlations that are important in high-T_c cuprates and oxypnictides. However, the most straightforward way to include correlation effects, the LDA + U method, fails in attempts to improve results of the DFT calculation for BaBiO₃. In the present work it is shown that the reason for this failure is a high oxidation state of the compound. A conventional basis set used in the LDA + U method is pure atomic orbitals, that is usually a good approximation for d-bands Wannier functions in transition metal compounds (with the important exception of high oxidation state materials). For BaBiO₃ one should use genuine Wannier functions as a basis set for the LDA + U correction potential, because they differ a lot from pure atomic Bi 6s-orbitals. Below it is demonstrated that the LDA + U correction in the Wannier function basis results in a strong improvement of the calculated parameters for crystal structure and spectral properties for BaBiO₃.

The LDA + U [1] and LDA + DMFT [2] methods for the electronic structure calculations of strongly correlated systems are designed as a combination of DFT and the model Hamiltonian approach (such as Hubbard or Anderson models). A DFT potential is defined by electronic density only and hence is invariant under unitary transformation of auxiliary Kohn–Sham orbitals (DFT does not consider the basis functions of a particular method; DFT considers only Kohn–Sham orbitals and total charge density). However, in model Hamiltonians, Coulomb interaction between electrons is explicitly defined for localized states having symmetry of atomic orbitals. A specific form of these orbitals remains to be determined in a calculation scheme developed for LDA + U method realization. For the LMTO [3] method, the basis set (Muffin-tin orbitals) is explicitly constructed in the form of atomic-like orbitals and so in this case the choice of orbitals for the Coulomb interaction term is obvious. For the LAPW [4] method, there are atomic symmetry wavefunctions inside muffin-tin spheres and there are pseudoatomic orbitals even for pseudopotential methods. These orbitals are usually used to define the Coulomb interaction term in a LDA + U calculation scheme based on corresponding methods.

The most general and 'rigorous' way to define atomic-like orbitals in a solid is to build Wannier functions. That is done via Fourier transformation from delocalized Bloch functions ψ_k in reciprocal space of wavevectors k to functions W^T in direct space lattice site vectors T. Both sets, the Bloch functions ψ_k and Wannier functions W^T, generate the same wavefunctions in Hilbert space and can be considered as different basis sets for this space. Correlation effects are usually strong enough to be taken into account for transition metal compounds, where d- or f-states produce narrow energy bands. Fourier transformation for the corresponding Bloch functions results in Wannier functions, which are well localized inside an atomic sphere and hence using atomic d- or f-orbitals instead of them is a good approximation.

Correlated states differ substantially from pure atomic d-orbitals in compounds with elements in a high oxidation state, for example, CaFeO₃,AuTe₂, and BaBiO₃. Narrow energy bands near the Fermi level have a serious admixture of p-states in such compounds, which could be even larger in value than d-states contribution. Moreover, an energy gap could be formed predominantly by p-states with a slight admixture of d-states (as it happens in AuTe₂) [5]. The mixed nature of correlated electronic states does not allow us to use traditional implementations of the LDA + U method with pure atomic orbitals as a basis set for the Coulomb interaction term. An explicit Wannier function basis set for correlated states should be chosen to perform electronic structure calculations for these materials.

Many researchers have used the Wannier function basis for strongly correlated materials calculations in recent years, for example, see [6–9]. In the present work we investigate BaBiO₃ where the Wannier function basis in contrast to atomic orbitals is crucial for proper crystal structure description. In the case of a BaBiO₃ partially filled band, which is a source of correlation effects, it is formed by Bi 6s-states and oxygen 2p-states. The energy of the original Bi 6s-level is below the oxygen band but strong hybridization between them results in the splitting of the partially filled band upward from the 2p-states manifold. The contribution of oxygen states to this band is stronger than the contents of Bi 6s-states but the corresponding wavefunction still has s-symmetry in respect to the Bi-site. The Wannier function constructed from this band has the same symmetry but it is spatially distributed more on oxygen sites than on the Bi ion. Accounting for the Coulomb interaction on such a Wannier function basis leads to results very different from the LDA + U calculations using pure Bi 6s-orbitals.

2.1. Wannier functions

Wannier functions (WFs) $\vert {W}_{n}^{\mathbf{T}}\rangle$ are defined as Fourier transforms of the Bloch functions |Ψ_nk〉:

$$\begin{eqnarray} &&\vert {W}_{n}^{\mathbf{T}}\rangle =\frac{1}{\sqrt{\Omega }}\mathop{\sum }\limits _{\mathbf{k}}{\mathrm{e}}^{-\mathrm{i}\mathbf{kT}}\vert {\Psi }_{n\mathbf{k}}\rangle , \end{eqnarray} \tag{ 1 }$$

where T is the lattice translation vector, n is the band number, and k is the wavevector. WFs seem to be an optimal choice to describe correlated states because these functions are localized, site-centered, and represent a complete basis set for the Bloch functions' Hilbert space. In the present paper WFs were generated as projections of Bloch sums of the atomic orbitals $\vert {\phi }_{n\mathbf{k}}\rangle ={\sum }_{\mathbf{T}}{\mathrm{e}}^{\mathrm{i}\mathbf{kT}}\vert {\phi }_{n}^{\mathbf{T}}\rangle$ onto a subspace of the Bloch functions (the detailed description of the WFs construction procedure within the pseudopotential method is provided in [10]):

$$\begin{eqnarray} \vert {W}_{n}^{\mathbf{T}}\rangle &=&\frac{1}{\sqrt{{N}_{\mathbf{k }}}}\mathop{\sum }\limits _{\mathbf{k}}\vert {W}_{n\mathbf{k}}\rangle {\mathrm{e}}^{-\mathrm{i}\mathbf{kT}}, \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} \vert {W}_{n\mathbf{k}}\rangle &\equiv &\sum _{i={N}_{1}}^{{N}_{2}}\vert {\Psi }_{i\mathbf{k}}\rangle \langle {\Psi }_{i\mathbf{k}}\vert {\phi }_{n\mathbf{k}}\rangle . \end{eqnarray} \tag{ 3 }$$

The resulting WFs have symmetry of the atomic orbitals ϕ_n and describe electronic states that form energy bands numbered from N₁ to N₂. To describe Coulomb repulsion between occupied and empty states one needs to compute occupation numbers of WFs. The WFs occupancy matrix ${Q}_{n m}^{\mathrm{WF}}$ is given by

$$\begin{eqnarray} &&{Q}_{n m}^{\mathrm{WF}}=\langle {W}_{n}^{0}\vert \left (\mathop{\sum }\limits _{\mathbf{k}}\sum _{i={N}_{1}}^{{N}_{2}}\vert {\Psi }_{i\mathbf{k}}\rangle \theta ({\varepsilon }_{i}(\mathbf{k})-{E}_{\mathrm{f}})\langle {\Psi }_{i\mathbf{k}}\vert \right )\vert {W}_{m}^{0}\rangle \end{eqnarray} \tag{ 4 }$$

where θ is the step function, ε_i(k) is the one-electron energy for the state i, and E_f is the Fermi energy.

2.2. LDA + U method in WF basis

One of the most general and accurate approaches to describe the electronic structure of strongly correlated systems is the dynamical mean-field theory (DMFT) [11–13]. The method combining the ab initio DFT approach and model DMFT calculations—LDA+DMFT was applied recently to the quantitative description of both the electronic and crystal structure of several compounds with strong electronic correlations [6, 14–17]. The LDA+DMFT method is successful and promising but it demands a lot of computational resources. In the present paper the static limit of dynamical mean-field theory, i.e. the LDA + U method, is used. There is an additional term in the Hamiltonian operator to take into account Coulomb correlations. If one considers on-site interaction only and takes into account the diagonal form of the occupation matrix ${Q}_{n m}^{\mathrm{WF}}$, the additional term is written in the WF basis as

$$\begin{eqnarray} &&{\hat {H}}^{\prime}(\mathbf{k})=\mathop{\sum }\limits _{m}\vert {W}_{m\mathbf{k}}\rangle (U\times (\frac{1}{2}-{n}_{m}))\langle {W}_{m\mathbf{k}}\vert , \end{eqnarray} \tag{ 5 }$$

where ${n}_{m}={Q}_{m m}^{\mathrm{WF}}$—the occupation number for the mth WF. In the case of BaBiO₃ the non-diagonal terms of ${Q}_{n m}^{\mathrm{WF}}$ are equal to zero by construction (we have strictly orthogonal WFs of Bi s-orbital symmetry). For more complex basis sets it could be necessary to use rotationally invariant potential correction and take into account non-diagonal terms.

The full Hamiltonian operator is written as

$$\begin{eqnarray} &&{\hat {H}}^{\mathrm{LDA}+\mathrm{U}}={\hat {H}}^{\mathrm{LDA}}+{\hat {H}}^{\prime}, \end{eqnarray} \tag{ 6 }$$

where ${\hat {H}}^{\mathrm{LDA}}$—the Hamiltonian operator within the LDA (or GGA) approximation. The form of the correction potential (5) results in a negative addition to potential  − U/2 for occupied states n_m = 1 and in a positive value  + U/2 for empty states n_m = 0. Hence the LDA + U correction to DFT will increase energy separation between bands below and above the Fermi level increasing energy gap values and enhancing the total energy gain of the distortion that causes the gap appearance.

The value of Coulomb interaction parameter U is computed via the constrained LDA calculation as described in [10].

The total energy is computed as

$$\begin{eqnarray} &&{E}^{\mathrm{tot}}={E}^{\mathrm{LDA}}+{E}^{U}-{E}^{\mathrm{DC}}, \end{eqnarray} \tag{ 7 }$$

where E^LDA is the total energy from a standard DFT calculation (LDA or GGA could be used, we preferred GGA), ${E}^{U}=\frac{1}{2}{\sum }_{m\not = m{^{\prime}}}U{n}_{m}{n}_{{m}^{\prime}},{E}^{\mathrm{DC}}=\frac{1}{2}U n(n-1)$ is the double counting correction, and $n=\mathrm{Tr}({Q}_{n m}^{\mathrm{WF}})$ is the total occupancy of WFs.

Undoped BaBiO₃ has a monoclinic crystal structure (symmetry group is C2/m). The structure could be obtained from an ideal cubic perovskite structure simultaneously breathing distortion of the BiO₆ octahedra and tilting the octahedra around the [110] axis. The primitive cell contains two formula units. Chemical formula could be written [18] as Ba₂²⁺Bi³⁺Bi⁵⁺O₆²⁻. Two different lengths of Bi–O bonds correspond to two different valences of Bi ion: the short to Bi⁵⁺–O and the long to Bi³⁺–O. The neighboring Bi-type ions for Bi³⁺ are Bi⁵⁺ and vice versa. The alternation of Bi³⁺ and Bi⁵⁺ ions forms a charge density wave. The experimental value [18] for the breathing distortion is b = 0.085 Å and for the tilting is t = 10.3°. Spectroscopy measurements show that BaBiO₃ is an insulator with the energy gap value  ≈ 0.5 eV  [19].

Liechtenstein et al performed calculations of the electronic and crystal structure of BaBiO₃ with the full-potential LMTO method [20]. The calculated values of the BiO₆ octahedra distortion (b = 0.055 Å,t = 8.5°) underestimate the experimental data as well as the energy gap value  ≈ 0.02 eV . Also authors were unable to obtain cubic crystal structure instability against pure breathing octahedra distortion without tilting.

Franchini et al have modeled the crystal structure distortions more successfully [21, 22]. The authors used the Heyd–Scuseria–Ernzerhof (HSE) hybrid functional that includes 25% of the exact Hartree–Fock exchange in addition to DFT, so it partially takes into account correlation effects. The hybrid functional allowed us to reproduce the energy gap value with a slight overestimation (0.65 eV) and get octahedra distortion values close to the experimental ones (b = 0.09 Å,t = 11.9°).

Thonhauser and Rabe [23] have reported the cubic crystal structure instability in respect of breathing distortion of the BiO₆ octahedra using the obsolete local density approximation (LDA) while the more rigorous generalized gradient approximation (GGA) calculations do not show any minimum at all (see figure 3). The total energy lowering value for distorted structure in LDA calculations [23] was 3.35 meV/f.u. corresponding to a crystal structure transformation temperature as low as 30 K, which clearly contradicts the experimental fact that the distorted cubic crystal structure of BaBiO₃ is stable at room temperature and higher temperatures [24].

In the present work to describe the correlated states with WFs, a self-consistent calculation within the GGA was performed as the first step. The pseudopotential plane-wave method implemented in the Quantum-ESPRESSO package [25], was used. Vanderbilt ultrasoft pseudopotentials [26] were taken from the QE pseudopotentials library. A kinetic-energy cutoff for the plane-wave expansion of the electronic states was set to 45 Ryd. Integrations in reciprocal space were performed using a (10, 10, 10) Monkhorst–Pack [27] k-point grid in the full Brillouin zone.

The GGA calculation results in a semiconducting solution with a very small indirect energy gap  ≈ 0.05 eV . Two energy bands in the interval [−1.5; 2.5] eV, see figure 2, are generally formed by the O p-states, but they have symmetry of the Bi 6s-states. The Bi 5d and Bi 6p partial densities of states give a vanishing contribution to the states near the Fermi level.

WFs were generated by projection of two Bi s atomic orbitals onto subspace defined by two energy bands near the Fermi level. Squared moduli of WF centered on the Bi³⁺ ion are shown in figure 1. WF has symmetry of the s-orbital. However, it has a significant contribution from p-states of neighboring oxygen ions. The second WF (centered on Bi⁵⁺ ion) will differ only by the ratio between O p and Bi s states.

Zoom In Zoom Out Reset image size

The contribution of oxygen states to partially filled bands could be accounted for by including in the Hamiltonian a Hubbard U correction to the oxygen p-states in addition to Bi s orbitals. However, in that case one needs to take into account also off-diagonal inter-atomic s–p terms in the occupation matrix. The WF approach proposed in the present work is simpler.

To perform calculation within the LDA + U method in WF basis, the U value for BaBiO₃ was computed via the constrained LDA method. The Hubbard U value could be calculated as the second derivative of total energy over occupation of the correlated state [28]:

$$\begin{eqnarray} &&U=\frac{{\partial }^{2}{E}_{\mathrm{DFT}}}{{\partial }^{2}{n}_{\mathrm{corr}}}. \end{eqnarray} \tag{ 8 }$$

In DFT the one-electron energy of a state is a derivative of total energy in respect to the state occupancy. Then the U value can be calculated as:

$$\begin{eqnarray} &&U=\frac{\partial {\epsilon }_{\mathrm{corr}}}{\partial {n}_{\mathrm{corr}}}, \end{eqnarray} \tag{ 9 }$$

where _corr is one particle energy of the correlated state, and n_corr—its occupation number. The partial derivative is taken numerically as described below. Constrained potential that affects only states considered as correlated ones is defined as a projection operator:

$$\begin{eqnarray} &&{\hat {H}}_{\mathrm{constr}}^{\mathbf{k}}=\mathop{\sum }\limits _{n}\vert {W}_{n\mathbf{k}}\rangle \delta {V}_{n}\langle {W}_{n\mathbf{k}}\vert , \end{eqnarray} \tag{ 10 }$$

where δV_n is a small number, for example, ±0.1 eV, positive for WF centered on one Bi-site and negative for WF centered on another Bi-site. The constrained potential is included in the Hamiltonian operator during the self-consistency cycle, which allows us to take into account screening of charge redistribution. As a result, the WF occupation numbers and energies change from unconstrained results. The derivative (9) is calculated as a ratio of these variations. We have checked that the calculated derivative value of U does not depend on δV_n, by using its various values equal to (0.1, 0.2, 0.3, 0.4 eV).

The resulting U value for BaBiO₃ equals 0.7 eV. The value of the Coulomb parameter strongly depends on the spatial form of the basis function [29]. Therefore the rather small value 0.7 eV obtained in our calculations is reasonable for WFs shown on figure 1 which is rather extended in space. The atomic Bi s-orbital is much more localized in space than WFs hence the U value for atomic orbitals was found to be significantly larger U = 2.5 eV .

The U value for the WFs basis is half the width of the partially filled energy bands, which indicates the relative weakness of Coulomb interaction in BaBiO₃. But, as is shown below, this Coulomb interaction plays an important role in the formation of the electronic states and crystal structure of the compound.

A simulation of the oxygen octahedra tilting and breathing distortion was performed in a monoclinic (I2/m) supercell with 2 f.u. Sixty primitive cells of Ba₂²⁺Bi³⁺Bi⁵⁺O₆²⁻ were generated corresponding to different combinations of the breathing and tilting distortions. For every structure a self-consistent calculation within GGA was done. Then WFs were constructed and self-consistent calculations within LDA + U (actually GGA+U) were performed. Coulomb interaction correction was applied only to the states described with the two WFs as defined in equation (5).

In figure 2 the formation of the energy gap due to the cubic structure distortion is shown. The band structure of ideal cubic BaBiO₃, see figure 2(a), is definitely metallic. Two energy bands in the [−1.5; 2.5] eV energy interval cross the Fermi level. These bands are double degenerate along the A → M direction. The cubic crystal structure is unstable against distortion that lowers the symmetry and opens an energy gap in band structure.

Zoom In Zoom Out Reset image size

In BaBiO₃ the breathing distortion of the oxygen octahedra plays the crucial role. In figures 2(b) and (c) band structure is shown for the crystal with the breathing-only distortions of the oxygen octahedra. The breathing value is b = 0.075 Å. The degeneracy breaking along the A → M direction is clearly seen in the GGA calculation result. Two separate bands appear in the [−1.5; 2.5] eV interval but the Fermi level crosses both of them and the energy gap does not open. Accordingly there is no total energy minimum due to the breathing distortion in the GGA calculation. That is clearly illustrated by the black curve in figure 3. The dependence of total energy of BaBiO₃ cell on the breathing distortion calculated within GGA has a minimum for the distortion parameter value b = 0 Å (stable cubic structure). Hubbard interaction correction (GGA + U) in the basis of WFs leads to the appearance of a band gap equal to 0.3 eV, figure 2(c). Also the cubic crystal structure of BaBiO₃ is unstable against the pure breathing distortion in the GGA + U calculation. The total energy of the cell has a minimum for b = 0.075 Å, see figure 3. The energy lowering is equal to 326 meV/f.u.

Zoom In Zoom Out Reset image size

The total energy dependence calculated with the GGA + U method in atomic orbitals basis is also shown in figure 3 with the dashed green curve. The minimum of the total energy corresponds to the cubic crystal structure. One can see that the atomic orbitals basis potential correction calculation does not result in an improvement of the crystal structure description. The energy dependence obtained with the LDA calculation is shown in figure 3 with the dashed magenta curve. It has a very shallow minimum for non-zero breathing distortion in agreement with [23], but there is no minimum in the curve obtained in GGA that is supposed to be a more accurate approximation than LDA.

In the monoclinic cell of BaBiO₃ with simultaneous breathing and tilting of BiO₆ octahedra the energy gap (equal to 0.05 eV) exists in the result of the GGA calculation, see figure 2(d). Despite the fact that the gap value is very small, the cubic structure is unstable against combined distortion in GGA. For the monoclinic structure the total energy minimum corresponds to the distortion parameter values (b = 0.04 Å,t = 12°). The values underestimate experimental data in agreement with the previous calculations [20].

The GGA + U correction for the monoclinic crystal cell increases the energy gap value to 0.55 eV, see figure 2(e), in good agreement with the experimental value of 0.5 eV.

The splitting of the half-filled double degenerate Bi s band and the opening of the energy gap results in a filled band corresponding to the 6s-state of one Bi-site (Bi³⁺) and an empty band formed by the states of another Bi-site (Bi⁵⁺). A charge difference between these two Bi-sites (equal to 2e) is expected for the fully ionic picture. However, the occupancy values difference for two atomic Bi 6s-orbitals in the LDA + U calculation is equal to 0.43e only in agreement with [21]. The small value of the charge difference can be understood from an analysis of WF contents. The Bi 6s-symmetry Wannier function for Bi³⁺ site has only 23% of Bi s atomic orbitals and 77% of O p orbitals; on Bi⁵⁺ site WF composition is 20% of Bi s and 80% of O p. A change in occupation number of WF by ±1 results in a change of s atomic orbital occupancy by +0.23 for Bi³⁺ and −0.20 for Bi⁵⁺, which gives 0.43 in total.

In figure 4 the total energy of BaBiO₃ in GGA and GGA + U as a function of breathing and tilting distortions is shown. The zero energy value corresponds to the GGA total energy of an ideal cubic crystal structure. The surface was obtained by extrapolation from 60 points as the fourth-order polynomial. The experimental values of the distortions are (b = 0.085 Å,t = 10.3°) and are marked by a cross in figure 4. The GGA calculations gave a total energy minimum for distortion value parameters (b = 0.04 Å,t = 12°) strongly underestimating the experimental data. The total energy minimum within the GGA + U calculation corresponds to distortion parameter values (b = 0.075 Å,t = 12°) that are very close to the experimental ones.

Zoom In Zoom Out Reset image size

Figure 4 clearly illustrates the effectiveness of the calculation scheme proposed in the present paper. The usage of the WFs basis for the partially filled states in the GGA + U method allows one to reproduce not only the energy gap value, but also noticeably improve the description of the crystal structure distortions in BaBiO₃. The most important result is the instability of the cubic crystal against the pure breathing distortion obtained in the GGA + U calculations.

The scheme for studying electronic correlations in high oxidation state compounds based on the LDA + U method in the Wannier functions basis was proposed and applied to the case of BaBiO₃. The energy gap value and crystal structure parameters of monoclinic BaBiO₃ were successfully described in agreement with experimental data. Instability of the cubic crystal structure in respect to the pure breathing distortion was obtained for the first time.

Support by the Russian Foundation for Basic Research under grants nos RFFI-10-02-00046a and RFFI-10-02-96011ural, the President of Russian Federation Fund of Support for Scientific Schools grant 1941.2008.2 and the Program of the Russian Academy of Science Presidium 'Quantum physics of condensed matter' is gratefully acknowledged. Dmitry Korotin acknowledges the support of the Dynasty Foundation.