Theoretical analysis of oxygen vacancies in layered sodium cobaltate, Na_xCoO_2−δ

Layered sodium cobaltate Na_xCoO₂ (x ≤ 1) has in recent years attracted considerable attention because of its remarkable electronic and magnetic properties. Most of the interesting physical properties of this material come from its quasi-two-dimensional CoO₂ layers where cobalt exists as Co³⁺ and Co⁴⁺ ions, surrounded by intercalating Na⁺ ions [1–3]. Understanding how the oxide layer is affected by temperature and sample composition is therefore crucial when interpreting the physical properties.

By changing the sodium content, one can continuously tune the oxidation state of the cobalt ions. In NaCoO₂ (x = 1) all cobalt atoms are found in a Co³⁺ state and the material is an insulator, whereas for x < 1, negatively charged sodium vacancies are introduced. Na_xCoO₂ is known to be a p-type conductor with metallic conduction behavior in which itinerant electron holes give rise to a particularly high thermoelectric effect [4, 5].

The influence of point defects other than Na vacancies on sodium cobaltate is still debated. Their concentrations are linked through the electroneutrality requirement and will thus vary depending on Na content, temperature and oxygen activity giving rise to a complex scenario. Among the possible defects, it has been recently reported that oxygen vacancies could have an interesting two-fold effect on the thermoelectric properties in this material. Firstly, they can affect the cobalt oxidation state, i.e. reducing the charge carrier concentration and increasing the thermopower. Unfortunately, such defects also modify phonon scattering so that the overall effect is a reduction of the figure of merit of Na_xCoO_2−δ with increasing the defect concentration (δ) [6].

Experimental data about the concentration of vacancies in sodium cobaltate are contradictory, especially at x ∼ 0.75. While some groups have identified an oxygen non-stoichiometry up to δ = 0.16 [6–11], others found no vacancies within the experimental error [12–14]. In order to clarify this, we have in this work studied the oxygen vacancy formation in Na_0.75CoO₂ by first principles density functional theory (DFT) calculations. From these results we computed free energies of formation, at various temperatures and oxygen partial pressures. We found that oxygen vacancies have high free energies of formation even at high temperature, suggesting that the experimentally observed weight losses attributed to oxygen vacancies may have a different origin.

We assumed the following reaction (using Kröger–Vink notation):

$$\begin{eqnarray} &&{\mathrm{Na}}_{x}{\mathrm{CoO}}_{2}\rightleftharpoons {\mathrm{Na}}_{x}{\mathrm{CoO}}_{2-\delta }+\delta {\mathrm{V}}_{\mathrm{o}}^{q \bullet }+q \delta {\mathrm{e}}^{\prime}+\frac{\delta }{2}{\mathrm{O}}_{2}. \end{eqnarray} \tag{ 1 }$$

The Gibbs' free energies of formation for oxygen vacancies ΔG_f(P,T) at a pressure P and temperature T were accordingly computed from [15]

$$\begin{eqnarray} \Delta {G}_{\mathrm{f}}(P,T)&=&{G}_{\mathrm{def}}(P,T)-{G}_{\mathrm{perf}}(P,T)\nonumber\\ &&+~\frac{1}{2}{\mu }_{{\mathrm{O}}_{2}}(P,T)+q{\mu }_{\mathrm{e}}. \end{eqnarray} \tag{ 2 }$$

In equation (2) G_def(P,T) and G_perf(P,T) are the Gibbs' free energy of the defective and pristine supercells respectively, q is the charge state per defect, and μ_i is the chemical potential of the species i.

Thermal contributions of the solid were considered as negligible (G_i(P,T) = H_i(P⁰,0)) with respect to those of the oxygen gas, whose chemical potential μ_O2(P,T) was calculated from

$$\begin{eqnarray} {\mu }_{{\mathrm{O}}_{2}}(P,T)&=&{H}_{{\mathrm{O}}_{2}}({P}^{0},0)+\Delta {\mu }_{{\mathrm{ O}}_{2}}({P}^{0},T)\nonumber\\ &&+~{k}_{\mathrm{ B}}T\ln \left (\frac{P}{{P}^{0}}\right ), \end{eqnarray} \tag{ 3 }$$

where k_B is the Boltzmann constant and Δμ_O2 is the change in the chemical potential when moving from T = 0 to T at constant pressure P⁰ = 1 bar, taken from thermodynamic tables. Enthalpies at standard pressure and T = 0 K,H_i(P⁰,0), were taken as the calculated DFT total energies and the corresponding concentrations of defects were then computed following:

$$\begin{eqnarray} &&\delta =n{\mathrm{e}}^{-\frac{\Delta {G}_{\mathrm{f}}(P,T)}{{k}_{\mathrm{B}}T}} \end{eqnarray} \tag{ 4 }$$

where n is the number of oxygen atoms per formula unit.

The system charge q was simulated by adjusting the total number of electrons in the supercell and at the same time adding a compensating jellium background to avoid diverging Coulomb contributions. The contribution of the background charge to the defect formation energy was taken into account by a correction term γ, here considered as the shift in the average electrostatic potentials at a bulk-like lattice site far from the vacancy in the defective (V_def) and pristine (V_perf) supercell [16]

$$\begin{eqnarray} &&\gamma ={V}_{\mathrm{def}}-{V}_{\mathrm{perf}}. \end{eqnarray} \tag{ 5 }$$

Then, this correction was included in the electron chemical potential, μ_e determined by the Na_xCoO₂ Fermi energy E_F [17].

$$\begin{eqnarray} &&{\mu }_{\mathrm{e}}={E}_{\mathrm{F}}+\gamma . \end{eqnarray} \tag{ 6 }$$

Periodic DFT as implemented in the VASP package [18] was used throughout this work. A spin polarized gradient corrected Perdew–Burke–Ernzerhof (PBE) functional [19] was used with a plane wave energy cutoff of 600 eV. Core electrons were included through the projector augmented wave (PAW) method [20]. Where needed we included an on-site Coulomb interaction (GGA + U) in order to localize cobalt d electrons [21]. The parameters for this system (U = 5.0 eV and J = 0.965 eV) were taken from [22] and are very similar to those used in other studies of the electronic structure of CoO₂-based oxides [1, 23, 24].

The reciprocal space was sampled by a Γ-centered k-point grid in which the maximum distance between points is 0.15 × 2π/|a|. Ionic positions were relaxed until the maximum force was lower than 0.03 eV Å⁻¹. Simultaneous relaxation of the lattice parameters was performed for the defect free structures.

To sample the many possible sodium ordered structures we adopted three different models for Na_0.75CoO₂. They are shown in figure 1 and named diamond, filled honeycomb (FHC) and zigzag. They are based on a tetragonal $\sqrt{3}\times 4\times 1$ supercell consisting of 16 formula units [25] and they represent different relative concentrations of the two inequivalent Na crystallographic sites, Na1 and Na2 (Na1/Na2 = 1,1/2, and 1/5). The three models are based on the experimental structure reported by Zandbergen et al [2] (diamond), and on the theoretical study of Meng et al [3, 26] (FHC and zigzag). These were chosen as a representative selection of simple periodic models from the large manifold of structures, suitable for a first principles study. Nevertheless, Na self-diffusion is an efficient process already at room temperature [27], so we expect a strong Na disorder at high T.

Zoom In Zoom Out Reset image size

The three supercells used as starting points for creating oxygen vacancies are shown in figure 1. The corresponding lattice parameters of our relaxed models (a = 11.45,b = 4.99,c = 10.88 Å) and the buckling of the CoO₂ layer (±0.05c) [28] are in good agreement with those determined by neutron diffraction experiments (a = 11.37,b = 4.92,c = 10.81 Å) [29]. Oxygen vacancies were considered in their neutral, +1 and +2 charge states (${\mathrm{V}}_{\mathrm{o}}^{\times },{\mathrm{V}}_{\mathrm{o}}^{\bullet }$ and ${\mathrm{V}}_{\mathrm{o}}^{\bullet \bullet }$ in Kröger–Vink notation), created at three different positions for each of the structural models. Oxygen atoms were removed from the sites shown as green balls in figure 1, then the whole lattice structure was relaxed to its equilibrium geometry. The electronic structure (density of states, DOS) of the perfect and defective Na_0.75CoO₂ (diamond structure) is shown in figure 2. Firstly we note that the electronic structures we computed are half-metallic, both for perfect and defective lattices, and independent of using GGA or GGA + U, in agreement with similar studies [23, 24]. The various sodium ion orderings and vacancy positions we have considered did not give any qualitative or quantitative differences in the DOS. Indeed Na states are found only at much more negative energies in the valence band, confirming that the intercalating ions simply provide charge carriers to the CoO₂ layers.

Zoom In Zoom Out Reset image size

For the pristine structure the dominant contributions to the DOS close to the Fermi level come from cobalt d and oxygen p states in the oxide layer, states that the distorted octahedral crystal field splits into ${\mathrm{t}}_{2\mathrm{g}}={\mathrm{e}}_{\mathrm{g}}^{\prime}\oplus {\mathrm{a}}_{1\mathrm{g}}$ manifolds [30]. Applying Bader charge analysis to our GGA results we found that all of the Co ions in the pristine structures are neither in a +3 or +4 state, but rather that holes are delocalized on all of the transition metal atoms. When the U term is considered, the holes are instead localized on few Co ions, whose d states are pushed to lower energies into the valence band, so that oxygen p states are now dominant at the Fermi energy (figure 2). However, the Hubbard-like U term is an arbitrary parameter used to correct GGA by applying an energy penalty to double occupation of d orbitals, for which many different values have been proposed ranging from 4.0 eV [31] to 5.5 eV [23]. While GGA + U was proven to localize holes in this class of materials it is not known how reliable the ground state energies (and their differences) are, in particular for metallic or half-metallic systems. Recently, Hinuma et al showed that GGA + U performs worse than GGA in reproducing experimental data for Na_xCoO₂ when x > 0.6, in particular its lattice parameters and the formation energy of different Na-ordered ground states [1]. Moreover, we note that within our approach the free energy of formation (equations (2) and (6)) is dependent on the position of the Fermi level of Na_0.75CoO₂. By introduction of the Hubbard-like term, the valence band edge (hence the Fermi level) shifts to lower energies, making the calculated free energy of formation strongly dependent on the arbitrarily chosen value of U. In GGA, the homogeneous oxidation state of Co ions is due to the self-interaction error that tends to over-delocalize the electron density. Still this scenario is compatible with recent NMR measurements which predict hole delocalization in Na rich cobaltates rather than localized Co⁺⁴ sites [32, 33]. It is also known that at low temperature different Na orderings establish patterns of Co³⁺/Co⁴⁺ which may modulate the formation energies of vacancies in neighboring oxygen sites, and this effect would not be correctly represented within our computational approach. Nevertheless, at high temperature the sodium sublattice is liquid-like [27] and fast diffusing Na ions would average out this effect. We then expect our GGA approach to be appropriate for representing qualitatively the high T phase, while it may fail representing the fine modulation induced by the spin states of cobalt ions at low temperature. Therefore, in this work we chose to rely only on free energies of formation computed with GGA, which are free from any arbitrary term.

By looking at the GGA electronic structure of the defective system in figure 2 we notice an enlargement of the band gap (of about 25% for the spin majority component), in which two sharp defect states lying at about 0.12 eV from each other are found. Charge density isosurfaces corresponding to the defect states are very similar in shape, and mostly localized on the vacancy nearest neighbors. Using a simple tight-binding argument for the bipartite CoO₂ lattice it is possible to predict how the removal of two p orbitals at the oxygen site, i.e. upon vacancy formation, as many degenerate states (known as midgap states) localize on the Co sublattice [34, 35]. In this case we suggest that the difference in symmetry and energy of the two defect levels may be caused by the Coulomb potential generated by the Na ions, which modulates the on-site energies of Co and O sites [36].

The analysis of the ground state electronic structure around the neutral defect ${\mathrm{V}}_{\mathrm{o}}^{\times }$ revealed that electron density accumulates at the vacancy site and on the three nearest neighboring Co ions. The difference in electron density for the neutral and doubly charged system is shown in figure 3. This difference in charge is too small to give a perfectly neutral defect site, but rather a partially positively charged vacancy surrounded by three Co ions which are reduced by 0.2 electrons each. We note here that GGA+U gives qualitatively the same results. The charge density at the vacancy site is removed by progressively increasing the system charge, resulting in an effectively doubly positively charged oxygen vacancy ${\mathrm{V}}_{\mathrm{o}}^{\bullet \bullet }$, as expected. The same analysis showed that Co neighbors have been reduced by 0.05e/atom when increasing the system charge to q = 2.

Zoom In Zoom Out Reset image size

Gibbs' free energies of formation (ΔG_f) for oxygen vacancies are shown in table 1 for the high temperature regime (T = 1000 K) relevant for applications as thermoelectric material. Results are reported for the three different models in both a $\sqrt{3}\times 4\times 1$ and $2\sqrt{3}\times 4\times 1$ supercell, corresponding to Na_0.75CoO_2−δ with δ = 0.031 and δ = 0.016, respectively. The table values were based on the Fermi energy E_F of the perfect material and P(O₂) = 1.0 atm.

Table 1.  Defect formation free energy as defined in equation (2) for neutral, singly and doubly charged oxygen vacancy in Na_0.75CoO_2−δ computed for two different δ values. T = 1000 K and P(O₂) = 1.0 atm. All the values are in eV.

	ΔGf (eV)
	δ = 0.031			δ = 0.016
	${\mathrm{V}}_{\mathrm{o}}^{\bullet \bullet }$	${\mathrm{V}}_{\mathrm{o}}^{\bullet }$	${\mathrm{V}}_{\mathrm{o}}^{\times }$	${\mathrm{V}}_{\mathrm{o}}^{\bullet \bullet }$	${\mathrm{V}}_{\mathrm{o}}^{\bullet }$	${\mathrm{V}}_{\mathrm{o}}^{\times }$
diam1	2.19	2.26	2.30	2.27	2.27	2.28
diam2	2.31	2.34	2.36	2.32	2.34	2.36
diam3	2.29	2.33	2.36	2.35	2.37	2.39
zigzag1	2.35	2.37	2.37	2.36	2.37	2.39
zigzag2	2.32	2.37	2.39	2.39	2.40	2.40
zigzag3	2.37	2.42	2.44	2.41	2.43	2.45
FHC1	2.46	2.44	2.48	2.40	2.39	2.42
FHC2	2.37	2.39	2.41	2.38	2.37	2.39
FHC3	2.19	2.24	2.26	2.20	2.23	2.25

Overall, the defect formation energies shown in table 1 are quite high, ranging from 2.19 to 2.48 eV at T = 1000 K and P(O₂) = 1.0 atm, implying a very low concentration of vacancies. Formation energies are distributed in a narrow range of 0.2 eV, with no clear preference for any of the structural models. This suggests that the Na ordering only plays a minor role for the vacancy formation energy. As well, the difference in the ΔG_f for neutral and charged defects is also small at both the concentrations δ considered, and lower than 0.1 eV in most of the cases. This suggests that the Coulomb potential generated by the charged defect is efficiently screened already at the length scale of the shortest lattice parameter, i.e. ∼5.0 Å. Nevertheless, it appears that the doubly charged defect ${\mathrm{V}}_{\mathrm{o}}^{\bullet \bullet }$ is slightly favored compared to the other two possible charges, suggesting a non-perfect dielectric screening. This small energy difference is a result of Na_0.75CoO_2−δ being a p-type conductor, with unfilled valence states. These states lie at energies only slightly above the Fermi level (see figure 2), and are readily filled by the extra electron(s) accompanying the singly charged and neutral oxygen vacancies.

In order to quantify the concentration of vacancies under relevant experimental conditions we calculated ΔG_f for x = 0.75 at different temperatures for the structure with the lowest free energy of formation (2.19 eV) in table 1. Results are shown in the left panel of figure 4 for three different temperatures (700, 1000 and 1300 K) as a function of P(O₂). In the right panel we show the corresponding concentration of vacancies per formula unit, δ. We also note that our formation energies are in good agreement with those recently reported by Yoshiya et al [37] for x = 0.50 and 1. It is clear that the oxygen vacancy free energy of formation is very high at relevant temperatures and oxygen partial pressures. Even at extreme experimental conditions such as T = 1300 K (close to the homogeneous melting point [38]) and P(O₂) = 10⁻⁷ atm the vacancy concentration would be as little as δ ≃ 6.5 × 10⁻⁴. If more common experimental conditions are chosen, e.g. T = 700 K and P(O₂) = 10⁻¹ atm, we obtain δ ∼ 10⁻¹⁸. This is several orders of magnitude lower than reported in some previous publications, e.g. [6].

Zoom In Zoom Out Reset image size

Having ruled out the Na ordering as the main influence on the vacancy formation energy we now test the role of sodium vacancies, thus indirectly that of the Co oxidation state. To do this we modify the sodium content x, as it is known to determine directly the magnetic behavior (hence the oxidation state) of Co ions by supplying electrons to the oxide layers [33]. We therefore studied the ${\mathrm{V}}_{\mathrm{o}}^{\bullet \bullet }$ defect formation for a range of Na concentrations, in order to span as much as possible the possible cobalt charges, from +4 (x → 0) to +3 (x → 1). Two different Na orderings were generated for each x value starting from the diamond structure, and the oxygen vacancy was created at three different sites in each of those. Averaged GGA results for all the structures considered are shown in figure 5, using the same T and P(O₂) conditions as in table 1. For Na rich structures, the change in the formation energies with x is substantial, but the absolute energies are still very large. When the Na content is reduced, ΔG_f lowers almost linearly down to x = 0.50, where it reaches a plateau of about 1.70 eV, suggesting that Co⁴⁺ ions favor the formation of O vacancies. However, the magnitude of the free energy of formation is still high, so that even at low sodium content (x ≤ 0.5) the concentration of oxygen vacancies under relevant experimental conditions is expected to be very low. A detailed analysis of the effect of oxygen vacancies onto the Co oxidation state or vice versa goes beyond the scope of this study.

Zoom In Zoom Out Reset image size

In conclusion, we have shown that the equilibrium concentration of oxygen vacancies in sodium cobaltate is very low at conditions relevant for applications as thermoelectric material. Our results may help to settle the disagreement between experimental reports on the oxygen stoichiometry in Na_xCoO₂. At high sodium content, concentrations of vacancies ranging from δ = 0.05 to 0.16 have been reported [6, 7, 10, 11], but other studies have concluded that the oxygen concentration is indeed stoichiometric [12–14]. Our results are clearly consistent with the latter ones.

We can here only speculate on the reason for this discrepancy. One possible explanation is existence of secondary phases in some of the experiments, where these phases are responsible for the oxygen weight loss upon reduction of oxygen partial pressure. Other possibilities comprise Na evaporation, the presence of other defects than considered here, or unintended interaction between the sample and experimental setup.

We acknowledge Tor Svendsen Bjørheim for helpful suggestions, the Research Council of Norway (Renergi project THERMEL-143386) for economic support, and the NOTUR consortium for providing access to their computational facilities.