Ab initio and empirical defect modeling of LaMnO_3±δ for solid oxide fuel cell cathodes

Abstract

Sr doped LaMnO₃ is a perovskite widely used for solid oxide fuel cell (SOFC) cathodes. Therefore, there is significant interest in its defect chemistry. However, due to coupling of defect reactions and inadequate constraints of the defect reaction equilibrium constants obtained from thermogravimetry analysis, large discrepancies (up to 4 eV) exist in the literature for defect energetics for Sr doped LaMnO₃. In this work we demonstrate how ab initio energetics and empirical modelling can be combined to develop a defect model for LaMnO₃. Defect formation enthalpies, including concentration dependence due to defect interactions, are extracted from ab initio energies calculated at various defect concentrations. Defect formation entropies for the defect reactions in LaMnO₃ involving O²⁻(solid) ↔ ½O₂(gas) + 2e⁻ are shown to be accessible through combining the gas phase thermodynamics and simple models for the solid phase vibrational contributions. This simple treatment introduces a useful constraint on fitting defect formation entropies. The predicted defect concentrations from the model show good agreement with experimental oxygen nonstoichiometry vs.P(O₂) for a wide range of temperatures (T = 873–1473 K), suggesting the effectiveness of the ab initio defect energetics in describing the defect chemistry of LaMnO₃. Further incorporating a temperature dependent charge disproportionation energy within 0.0–0.2 eV, the model is capable of describing both defect chemistry and oxygen tracer diffusivity of LaMnO₃. The model suggests an important role for defect interactions which are typically excluded from LaMnO₃ defect models, and sensitivity of the oxygen defect concentration to the charge disproportionation energy in the high P(O₂) region. Similar approaches to those used here can be used to model the defect chemistry for other complex oxides.

---

A. Introduction

Sr-doped LaMnO₃ (LSM) is presently used as a cathode materials in commercial zirconia electrolyte based solid oxide fuel cells (SOFC) due to its good electrical conductivity, good stability, reasonable catalytic activity for oxygen reduction reaction (ORR), and relatively low cost^1,2 Defects in perovskites have been shown to play a crucial role in ORR. For example, Ba_1−xSr_xCo_1−yFe_yO₃, one of the most active perovskites for ORR, has oxygen vacancy concentrations of near 23% under SOFC conditions.³ More broadly, many aspects of oxygen kinetics and catalysis on LSM are intimately related to its defect chemistry. Oxygen transport is governed by the oxygen tracer diffusion coefficient (D*), which is approximately proportional to oxygen vacancy content. In addition, the log of the surface exchange coefficient, k*, which governs how quickly oxygen can enter the material, is typically proportional to log of D* with a slope of approximately 0.5 (i.e., k* is proportional to the square root of D*) in mixed electronic-ionic conducting (MEIC) perovskites,^4–6 including LSM. The correlations between k* and D* implies oxygen vacancies play an important role in determining the magnitude of k*. Finally, the ORR activity is strongly dependent on k* and D* values,⁷ and hence also critically affected by oxygen vacancy content. These observations indicate the importance of understanding defect chemistry of perovskite as SOFC cathode catalysts.

Experimental efforts to directly obtain defect concentrations face some significant challenges. In particular, due to low oxygen vacancy content in LSM under typical SOFC conditions, direct measurement of bulk point defect concentration remains very difficult. Thermogravimetry (TG) analysis can be used to determine overall oxygen non-stoichiometry (δ in (La,Sr)MnO_3±δ), although such analysis does not independently identify the extent of cation vs.oxygen vacancies. To obtain detailed relationships of defect concentration vs.P(O₂) and temperature, defect models^8–16 are generally used, with experimental data such as oxygen nonstoichiometry as constraints. Such defect models can be used to solve for the concentrations of defect species in LSMvs.P(O₂) and temperature.

However, due to uncertainties in the TG determined equilibrium constants (up to 1–2 order magnitude, as shown in ref. 11) as well as the associated uncertainties in defect reaction enthalpies determined from the equilibrium constants (up to eV scale), significant discrepancies exist in the defect energetics among the previously developed LSM/LaMnO_3±δ (LMO) defect models. Defect reaction enthalpies can vary by 2–4 eV between different models, as summarized in Table 2 below. The lack of robust constraints from the TG data in the LSM/LMO defect models leads to multiple reaction energies that all show good agreement with the experiments, but leave significant uncertainty about the true values of these key defect energetics.

Non-ideal behavior is a possible source of significant uncertainty is LSM/LMO defect model parameters. Since defect interactions (non-ideality) are found to be substantial in LMO/LSM,^13,15 the constant defect reaction enthalpies and entropies derived from Van’t Hoff plots (equilibrium constants vs. 1/T plots), which neglect non-ideality of defects, will lead to significant errors at lower temperatures. The assumption of ideality may be the source of the disagreement generally found in the LMO/LSM random defect (ideal solution) models when utilizing defect reaction enthalpies and entropies obtained at high temperature to fit with the low temperature TG data, especially for the oxygen excess region. To include non-ideality in LMO and LSM defect chemistry, Mizusaki et al.¹³ proposed a modified defect model based on the observation of maximum oxygen overstoichiometry of LMO/LSM at high P(O₂) and applied an additional constraint to account for the saturation of cation defects, as explained based on the concept of vacancy excluding space around cation vacancies and Sr dopants.¹³ A consequence of this simplified interaction model is that the Mizusaki model predicts an abrupt flattening of O overstoichiometry for LMO/LSM at low temperatures as the cation vacancy concentration reaches the saturation limit, which implies a step function for the cation defect interaction from 0 to ∞ when the cation defects become saturated. Recently, Mebaneet al.¹⁵ modified the random defect model with a defect concentration dependent strain parameter (to account for the excess free energy from the strain field interactions between cation defects) and fit with experimental O nonstoichiometry of oxygen excess La_0.9Sr_0.1MnO_3±δ and La_0.8Sr_0.2MnO_3±δ over a wide range of temperatures (T = 873 K–1273 K). Despite the physical motivation and extensive fitting in Mebaneet al.'s strain defect model, it is not clear how well constrained the parameters are or if other interactions could play a role.

In this work, we seek for guidance from first principles calculations to resolve uncertainty in the defect energetics of the LMO/LSM defect models. An ab initio and empirical defect model is developed by explicitly calculating the vibrational entropies of the defect reactions as well as incorporating the calculated ab initio defect energetics, to serve as new constraints for the defect model. Ideally an ab initio defect model would use ab initio methods to calculate energies and entropies for all possible defect reactions to establish the defect model. However, the considerable number of defect species along with all possible coupling/interactions between the defects in multicomponent oxide systems such as perovskites generally produce enormous amount complexity. In the case of LMO/LSM, the number of the defect species (vacancies, antisite defects, and different valence states of Mn) is more than 10 and the number of the defect coupling terms is around 100, hence including all possible defect reactions and interactions in the defect model is not practical. Furthermore, validation of such a complex model would be very difficult if not impossible with present experimental data. Therefore, in this study we firstly justify the assumptions made in most of the previous defect models^8–16 by calculating reaction energies of possible defect reactions and identifying the important defect reactions that should be included in the LMO defect model. Once the key reactions and interactions have been identified, ab initio methods have been utilized to propose new values for the key reaction energies. We hope that the success of the present model will motivate further work to expand the number of reactions and defect interactions in complex oxides treated with fully ab initio based approaches.

Focusing on LMO as an example system, this work serves as a first step to demonstrate that a simple defect model and calculated ab initio defect energetics can predict defect concentrations which agree well with the experimental LMO nonstoichiometry over a wide range of temperature (873–1473 K). Overall, the model developed in this work allows one to directly utilizing ab initio defect energetics to improve understanding of defect chemistry of electron-rich perovskites as SOFC cathode materials. The derivation of the defect reaction vibrational entropies that involve exchange of oxygen ions with the oxygen gas phase can provide quantitative understanding of the vibrational free energy contributions at high temperatures for oxide native defect formation reactions. We readily accept that at this stage ab initio defect energies, particularly for complex oxides, come with their own significant limitations in accuracy. However, the approach demonstrated here is remarkably successful for LMO and can be applied to other systems in combination with experiments and ever improving ab initio methods for more accurate models.

This paper is organized as the following. Section B discusses which defect reactions are important and should be included in the defect model based on the calculated ab initio defect reaction energies. Section C describes the expressions of the LMO defect model with derivation of the vibrational entropy/free energy for the defect reactions. Section D discusses the ab initio modeling approach and compares the calculated LMO defect energetics with the reported defect energetics in the literature. The predicted oxygen nonstoichiometry, defect concentrations vs.P(O₂) at various temperatures, and oxygen tracer diffusivity, are discussed in Section E. Finally, conclusions are summarized in Section F.

B. Justification of important defect reactions

To justify which defect reactions are important and should be included in the LMO defect model, we firstly calculate 12 defect reaction energies based on total energies of perfect and defect containing 2 × 2 × 2 (40 atom) perovskite supercells and La and Mn binary oxides (the details of ab initio calculation approach are described in Section D below). To account for the temperature and oxygen partial pressure effects, we further combined the calculated ab initio defect reaction energies with the vibrational entropy contribution (as described in Section B.2 below) and an oxygen partial pressure term, ½·k_b·T·ln(P(O₂)), for the defect reactions involving O²⁻(solid) ↔ ½·O₂(gas) + 2·e⁻, to derive the corresponding reaction free energies without the configurational entropy contribution at T = 1000 K. For the rest of the defect reactions that do not involve oxygen exchange between the solids and the O₂ gas phase, the calculated ab initio reaction energies are used directly to approximate the . The results of the vs.P(O₂) at T = 1000 K for the 12 defect reactions are shown in Fig. 1. These defect reactions are given in Kröger-Vink notation and are based on small polaron model expressions. Here A and B refer to the La and Mn sublattices, respectively. It is noted the all the calculated do not include nonideality, but the ideal solution analysis allows identification of important defect reactions to be further used in the defect model.

Fig. 1 The calculated 12 defect reaction free energies without configurational entropy contribution vs.P(O₂) at T = 1000 K. The three major reactions to be further included in the defect model are (the three solid lines): Reaction (1) (oxygen vacancy formation reaction), reaction (2) (cation vacancy formation reaction), and reaction (12) (charge disproportionation reaction). Reaction (2) is used to represent reactions (3)–(7). Reaction (9) is also neglected as the estimated Mn antisite defect concentration based on the of reaction (9) is small (∼10⁻⁴).

Based on the obtained , the defect reactions which are important to be included in the defect model are: reaction (1) (oxygen vacancy formation), reaction (2)–(7) (cation vacancy formation), reaction (9) (Mn antisite defect formation on La sublattice), and reaction (12) (charge disproportionation). These reactions are identified as important since these are significantly lower than the others for at least some PO₂ values. Among the various cation vacancy formation reactions (reactions (2)–(7)), the B site cation vacancy formation reaction is about 0.2 eV more stable than the A-site cation vacancy formation reaction in the high P(O₂) region, which is where cation vacancies can form in significant concentrations. This 0.2 eV energy difference corresponds to about a factor of 10 difference in concentrations at T = 1000 K, which suggest B-site cation vacancies are energetically more favorable to form than the A site vacancies. However, in the defect model below (see Section C) we consider the use of a combined A-site and B-site cation vacancy defect formation reaction (i.e.reaction (2)) to represent the various cation vacancy formation reactions. There are a few reasons to make this approximation. First, this treatment leads to the constraint of equal concentration on the A and B sublattices, which provides a significant simplification in the modeling. Furthermore, the small differences in the cation vacancy formation energies which we are ignoring will not cause significant change in the predicted overall cation vacancy concentrations at high P(O₂), the oxygen nonstoichiometry, or the oxygen diffusivity. Finally, it is not clear that including the different cation occupancies with our present ab initio numbers would actually lead to a more accurate model. The of reactions (2)–(7) are so close that it is difficult to be sure the distinction is meaningful as the uncertainties in the energies used in this work are at least ∼0.1 eV. In addition, the experiments are neither consistent with the modeling results nor with each other. The A-site vs. B-site vacancy concentration ratio measured by neutron diffraction reported in ref. 17 is about 1.8, which suggests the defect energetic difference between the A-site and B-site cation vacancy formation is about 0.05 eV at 1000 K. However, other studies have found equal A-site and B-site vacancy,¹⁸ and up to 10% A-site vacancy vs. filled B-site.¹⁹ This spread in results makes it hard to validate the calculations against experiments. It is clear that additional computational and experimental work is needed to help resolving the cation vacancy formation energetic difference and we therefore use a combined A-site and B-site cation vacancy defect formation reaction.

The one other defect reaction of low enough energy to potentially play a role is the Mn antisite defect formation (reaction (8)). While this reaction has a low energy it requires a cation vacancy to be present on the A site, and therefore contributes only a small modification to that already small cation vacancy concentration. The estimated Mn antisite defect concentration based on the of reaction (8) is about ∼10⁻⁴ at 1000 K, which is very small and can be neglected.

Overall, these ab initio calculations suggest three major defect reactions that should be included in the LMO defect model: reaction (1): oxygen vacancy formation reaction, reaction (2): cation vacancy formation reaction, and reaction (12): charge disproportionation reaction. These reactions have been used previously and the ab initio results therefore validate the previous defect model assumptions.^8–16

C. Defect modeling approach

The following sections describe the details of the LMO defect modeling approach based on the calculated ab initio defect energetics:

1. The three major LMO defect reactions

As discussed above in Section B, the three independent defect formation reactions in the LMO defect model of this work are:

(i) Oxygen vacancy formation reaction (reduction reaction):

(ii) Cation vacancy formation reaction (oxidation reaction): (iii) Charge disproportionation reaction:It is noted that the Schottky defect reaction can be derived from linear combination of the above reactions , i.e., 3 × eqn (1) + eqn (2) − 6 × eqn (3) (thus is not a independent defect reaction in the defect model).

2. Expression of equilibrium constants for the assumed LMO defect reactions

For an oxide defect reaction at equilibrium, the reaction equilibrium constant can be written as:where n_j is a stoichiometry coefficient of a defect species j, is a chemical potential without the ideal configurational entropy of a defect j, k_b is the Boltzmann constant, T is temperature, and c_j is concentration of a defect j. is defined as free energy of a defect reaction without configurational entropy at P(O₂) = 1 bar and is a function of defect concentration (due to defect interactions) and temperature (due primarily to vibrational free energy contributions). Note that while the enthalpy term in will be allowed to have non-ideal terms, the entropy will be assumed to be the ideal entropy. This approximation has been widely used in CALPHAD approaches and yielded excellent accuracy for many systems.^20,21

3. Derivation of vibrational free energy/entropy for the LMO defect reactions

To derive an expression for of the LMO defect reactions (eqn (1)–(3)), we split into two parts: the defect reaction energy as a function of defect concentration, Δε_defect form(c_defect), and the vibrational free energy contribution of the defect reactions, ΔG^vib_{defect form}(T). As we are aiming to have defect energetic accuracy at ∼0.1 eV scale at ∼1000 K, the vibrational degrees of freedom are approximated to be the only non-configurational degrees of freedom of the defect reactions (other non-configurational degrees of freedom contributions, e.g., magnetic, are neglected under the assumption that their contributions tothe free energy are either small or cancel in the reactions). For a defect reaction involving O²⁻(solid) ↔ ½·O₂(gas) + 2·e⁻, we assume the vibrational free energy of the solid phase atoms not involved in the reaction do not change significantly during the reaction. Consequently, the main contribution of vibrational free energy for the defect reactions will be from the vibrational free energy difference between the O₂ gas and the exchanged O²⁻ ions of the solid. With this assumption one can calculate ΔG^vib_{defect form}(T) by incorporating the empirical thermodynamic free energy/entropy data²² for the O₂ gas phase and a simple Einstein model for the exchanged O²⁻ ions in the LMO solid phase:^23–26here x is the number of O involved in the defect reactions (x = 1/2 for the oxygen vacancy formation reaction, i.e.eqn (1), and x = −3/2 for the cation vacancy formation reaction, i.e.eqn (3)), and are the enthalpy and entropy, respectively, of the O₂ gas phase at P(O₂) = 1 bar (per O₂) (taken from the NIST thermodynamic database²²), and is the empirical free energy of the oxygen gas phase at P(O₂) = 1 bar. G_vib(O²⁻_solid) is the vibrational free energy for an oxygen ion in LMO and H⁰_vib(O²⁻_solid) is the vibrational enthalpy of an oxygen ion in LMO at T = 298.15 K. The subtraction of H⁰_vib(O²⁻_solid) from G_vib(O²⁻_solid)is to set H⁰_vib(O²⁻_solid) at 298.15 K as the reference ( is also referenced to ). Following ref. 25, G_vib(O²⁻_solid) and H⁰_vib(O²⁻_solid) can be calculated by using a simple Einstein model with the Einstein temperature (θ_E) equal to 500 K for oxygen ions in the solid phase LMO. The value θ_E = 500 K was determined by diagonalizing the local force constant matrix for O (all other atoms fixed) in LMO from ab initio calculations. This value is consistent with experimental Raman frequencies, which range from 140–660 cm⁻¹ (which corresponds to θ_E = 202 K–880 K).²⁸ The second term of the right hand side eqn (5) allows one to calculate the vibrational free energies of LMO defect reactions involving exchange of oxygen between LMO solid and O₂ gas by using the expressions shown in Table 1 for each term in the bracket of the eqn (5).

Table 1 Expressions of the terms in the bracket of eqn (5) for calculating vibrational free energy of a defect reaction involving exchange of oxygens between the solid and O₂ gas phase^23–26

Expression

a t = T/1000, A = 3.113 × 10⁻¹, B = 9.092 × 10⁻², C = −4.133 × 10⁻², D = 8.170 × 10⁻³, E = −7.686 × 10⁻³, F = −1.174 × 10⁻¹, G = 2.448 from ref. 22. b k _b = 8.617 × 10⁻⁵ eV K⁻¹, θ_E = 500 K,²⁵ and T_R = 298.15 K from ref. 25 and 27.

^{22 a}

^{22 a}

^{25,27 b}

^{25,27 b}

---

To demonstrate the accuracy of this simple model, the model entropy terms and experimental fits to entropies of reactions have been compared. Note that we do not discuss the predicted enthalpy terms. Within a harmonic model these will be set by the equipartition theorem²⁹ to approximately 3k_bT at the temperatures of interest. If we assume that anharmonic effects make a relatively small contribution then the proposed Einstein model will capture the enthalpy contributions accurately. In Fig. 2, we show our calculated vibrational entropies between T = 800–1300 K for O₂ gas (, the dotted a line), oxygen ions in the solids (, the dashed lines), and any defect reaction involving exchange of oxygens between O₂ gas and solid (Δs, the solid lines). It is seen in Fig. 2 that the dominant vibrational entropy term is from the O₂ gas phase, while the vibrational entropy for an O²⁻ in LMO is about 30% of the O₂ gas entropy. Nonetheless, at SOFC temperatures (∼1000 K), the vibrational entropy for oxygen ions in LMO will lead to ∼0.4 eV/O contribution, which will cause significant error in the calculation if not included. Overall, the calculated defect reaction entropy is about 80 J/(mol K) per O²⁻(solid) ↔ ½·O₂(gas) + 2·e⁻ between T = 800–1300 K (positive for the forward reaction and negative for the backward reaction), which is in good agreement with the reported reaction entropies of LMO defect reactions (78 J/(mol K) per O in ref. 8 and 88.7 to 96 J/(mol K) per O in ref. 10). Fig. 2 also shows that a variation of 100 K in the estimated θ_E for O²⁻_solid only causes a few J/(mol K) shift in the defect reaction entropy, corresponding to a change of reaction free energy in the range of 0.01–0.1 eV at T = 1000 K. This result suggests an error in the adopted Einstein characteristic temperature will not lead to significant change in defect reaction free energies and justifies the use of such a simple model for the solid phase oxygen vibrational thermodynamics.

Fig. 2 Calculated vibrational entropy terms for the O₂ gas phase (; the dotted a line), oxygen ions in the solids (; the dashed lines), and any defect reactions involving exchange of oxygens between O₂ gas and solid (Δs, the solid lines) based on the expressions shown in Table 1. The thick solid and dashed lines are for the results with θ_E = 500 K and the thin solid and dashed lines are for the results of θ_E = 400 K (in red) and 600 K (in green).

As the charge disproportional reaction does not involve exchange of species from the solid to gas phase, the magnitude of the reaction vibrational free energy and entropy is expected to be much smaller than the other two defect reactions discussed above. Therefore the vibrational entropy change during the charge disproportionation reaction is assumed to be zero for LMO and only the configurational entropy and the reaction energy are included.

It is also worth pointing out that the reaction vibrational entropy of the Schottky defect reaction will be equal to zero following the above approach, since the overall reaction does not involve exchange of oxygens between the solid LMO and the O₂ gas phase.

4. Expressions of LaMnO₃ defect energetics

With the explicit treatment for the vibrational free energy and entropy of the LMO defect reactions described in Section C.3, the only unknowns left in the LMO defect model are defect reaction energies (Δε_defect form(c_defect) in eqn (5)). To include non-ideality of LMO defect reactions with respect to a range of defect concentration (oxygen and cation vacancies), the defect reaction energetics can be expressed as polynomials in defect concentration. Here, we include defect interactions between cation vacancies and interactions between oxygen vacancies but no interaction between the cation and oxygen vacancies. Although oxygen and cation vacancies may have interactions they are not present together in significant quantities at any P(O₂) or temperature in LMO (see Fig. 7) and their interactions can therefore be ignored.

In the following, we perform a series of first-principles calculations to obtain the LMO defect energetics. We use first order polynomials of defect concentration to represent the defect reaction energetics for the oxygen and cation vacancy formation reactions. A constant charge disproportionation reaction energy is also derived from the ab initio modeling. These energies will be discussed in Section D below.

D. Ab initio defect energetics

1. Ab initio modeling approach

Calculations were preformed with the Vienna Ab initio Simulation Package (VASP)^30,31 using Density Functional Theory (DFT) and the Projector-Augmented plane-Wave (PAW) method.^32,33 Exchange–correlation was treated in the Perdew-Wang-91 (PW-91)³⁴ Generalized Gradient Approximation (GGA) with electronic configurations of La (5s2 5p6 6s2 5d1), Mn_pv (3p6 3d6 4s1), and O_s (soft oxygen pseudopotential, 2s2 2p4). The energy cutoff was set to 450 eV and the brillouin zone was sampled by a Monkhorst–Pack k-point mesh of 8 × 8× 8 for a primitive ideal perovskite unit cell. Energy convergence with respect to k-points was better than 3 meV per perovskite formula unit. Bulk calculations are simulated with 40-atom (2a_c × 2a_c × 2a_c; a_c is the lattice constant of an ideal cubic LaMnO₃ perovskite), 80-atom (2a_c × 2a_c × 4a_c and ), 160-atom (2a_c × 4a_c × 4a_c and ), and 320-atom (4a_c × 4a_c × 4a_c) supercells with ions relaxed internally to restore BO₆ octahedron distortion. The above supercells implicitly neglect the Jahn–Teller distortions.²⁵ However, these are negligible at the high working temperature appropriate to SOFC (T = 800–1300 K) as the Jahn–Teller transition temperature, T_JT, is about 750 K. The LMO will therefore exhibit bulk structures with an averaged cubic symmetry under SOFC conditions. For the simulated magnetic structure, we use the ferromagnetic (FM) state to approximate the conductive characteristics of LMO at SOFC operating temperatures (above T_JT).

To treat the electron correlation errors in LMO using DFT with GGA, we adopt the GGA + U method,^35,36 where U_eff (=Coulomb U − exchange J) is applied to the Mn d orbitals, to improve accuracy of redox reaction energetics in our calculations.^25,37 Since it has been shown the redoxU_eff for Mn binary oxides are slightly different for different Mn oxidation states,³⁷ we adopt U_eff = 4 eV for calculating energies of the oxygen vacancy containing LMO (reduction of Mn³⁺ to Mn²⁺) following ref. 25 and 37, and U_eff = 3.5 eV for calculating the cation vacancy containing LMO (oxidation of Mn³⁺ to Mn⁴⁺). The latter value is chosen by fitting to the experimental oxidation enthalpy of 2Mn₂O₃ + O₂ → 4MnO₂ (1.678 ± 0.075 eV/O₂³⁸) with the approach suggested in Ref. 37. The different U_eff for Mn oxidation and reduction is also applied in the charge disproportionation energy (E_{chg_disp}) calculations, which are derived via the following equation:

E_{chg_disp} = [E_LMO(U_eff = 3.5 eV) + E_LMO(U_eff = 4 eV)] − [E_LMO+h˙(U_eff = 3.5 eV) + E_LMO+e′(U_eff = 4 eV)] (6)

where E_LMO is the calculated total energy of the perfect LMO, E_LMO+h˙ and E_LMO+e′ are the calculated total energies of stoichiometric LMO bulk with addition of one background charge in the simulated supercells, and the applied U_eff values are expressed in the parentheses.

As both the oxygen vacancy formation and the cation vacancy formation reactions in LMO are involved with exchange of oxygen between the solid phase and the O₂ gas phase, an O₂ gas energy reference is needed for calculating the defect reaction energies. In DFT ab initio modelling, the O₂ energy reference can be calculated by combining the ab initio energy of an isolated O₂ molecule along with a correction term to bring the experimental data and the quantum-mechanical calculations into agreement.^23,25,37,39 The applied O₂ energy correction term (eV/O₂) used here is 0.33 eV/O₂²⁵ and was determined for the exchange correlation functional (PW-91) and oxygen pseudopotential (soft oxygen pseudopotential) used in this work. In the following section, the combined energy, , which is −8.76 eV/O₂ as reported in ref. 25, will be used as the O₂ energy reference.

2. Ab initio LMO defect reaction energies vs. defect concentrations

In this work the relationships of defect formation enthalpies vs. defect concentration are obtained from supercells with different sizes. It is assumed that for a given defect concentration the energy is only weakly dependent on defect configuration in the relevant limit of low defect concentration (approximately δ = 0–0.2 for the oxygen excess LaMnO₃₊δ and δ = 0–0.125 for oxygen deficient LaMnO_3−δ). This assumption was checked by calculating different defect arrangements at fixed composition through use of different defect supercells. It was found that the calculated total energies of the simulated defect containing supercells (with either a separated A-site and B-site cation vacancy pair or an anion vacancy) with different supercell shape (2a_c × 2a_c × 4a_c and for 80-atom supercells, and 2a_c × 4a_c × 4a_c and for the 160-atom supercells) are fairly close in energy (within 0.1 eV per supercell). This result suggests that, provided the defects remain relatively far apart, the distances between defects are large enough to minimize the configurational dependence of the defect interactions in the calculations presented here.

It has been reported that vacancies in transition metal perovskites can cause chemical expansion^40,41 and the changes in lattice parameters with defect concentration are consequently of some importance. Therefore, in Fig. 3, we show lattice constant percent strain (ε = a_c(δ)/a_c(δ = 0)) vs. the corresponding defect concentration for the oxygen and cation vacancy containing supercells. Volume changes from the smallest defect concentrations (320 atom cells) yield formation volumes of −15.37 ų and 6.66 ų for a single cation (averaged over A and B sites) and anion defect, respectively. It is interesting to note that the cation vacancies cause lattice contraction, opposite the trend observed for the anion vacancies.

Fig. 3 Lattice constant percent strain of the calculated supercells at various defect concentrations: red solid diamonds (the solid line) represent the lattice constants of oxygen vacancy containing supercells, while blue solid circles (the dashed line) represent the lattice constant of the cation vacancy containing (1 A-site and 1 B-site vacancies) supercells. The calculated lattice constant of the perfect bulk LMO in this work is 3.97 Å.

Fig. 4 shows the calculated ab initio total energies of LMO (per Mn) referenced to that of the perfect bulk at the sampled O vacancy and cation vacancy concentrations. Note that for cation vacancies an A and B site vacancy are always created together. For simulating the cation vacancy formation reaction (eqn (2)), an A-site and B-site cation vacancy pair is modelled by putting the A-site vacancy at the corner of the supercell and B-site vacancy is placed at the farthest site from the corner. The 2 × 2 × 2 supercell data for cation vacancy formation reaction is not included in Fig. 4 because the simulated cation vacancy concentration (12.5% for both A and B site vacancy) is much higher than the maximum of experimental values (∼5%) at SOFC conditions. Furthermore, at such high cation vacancy concentration, the 2 × 2 × 2 supercell with 1 A-site and 1 B-site vacancies always has a first nearest neighbor cation vacancy pair, which could have strong short range configurational defect interaction. In such case, configurational dependency needs to be taken into account.

Fig. 4 Quadratic fittings of ab initio total energies (normalized as per Mn and referenced to the perfect bulk) vs. defect concentration calculated with various supercell sizes for (a) oxygen vacancy formation reaction (i.e.eqn (1)) and (b) cation vacancy formation reaction (i.e.eqn (2)). These allow to extracting relationships between defect reaction energies as a function of defect concentration.

In Fig. 4, the calculated total energies of the defect containing LMO at various defect concentrations, ε^total_{per Mn}(c_defect) (normalized as per Mn), are fitted with 2nd order (quadratic) equations of defect concentration:

ε^total_{per Mn}(c_defect) = ε₀ + ε₁·c_defect + ε₂·c²_defect (7)

where ε₀, ε₁,and ε₂ are zeroth, first and second order fitting parameters (note that ε₀ is equal to 0 by setting the energy of the perfect LMO bulk as the reference). The defect reaction formation energy (ε_defect form in eqn (5)) associated with a discrete change in defect concentration, δc_defect, can be written as:where x is defined as in eqn (5), and ε_O2 is discussed in Section C.1. In the last step we take the limit of δc_defect approaching zero to obtain the defect formation energy.

Rearranging eqn (8), we obtain the equation of defect reaction energy vs. defect concentration as:

where E_rxn = ε₁ + x·ε_O2/2, and I_rxn = 2·ε₂. The extracted E_rxn's and I_rxn's from Fig. 4 are further summarized below in Table 2, where E_red and I_red are for the oxygen vacancy formation reaction and E_ox and I_ox are for the cation vacancy formation reaction. In Section D.4 below, these defect energetic parameters (E_rxn's and I_rxn's) are further discussed and compared with those of the previous LMO defect models in literature.

Table 2 Comparison of the ab initioLaMnO₃ defect reaction energetics in this work and the reported values in literature

Defect reactions	Defect reaction energies vs. cdefect Δεrxn(cdefect) = Erxn + Irxn·cdefect, (eV)
	Ab initio (this work)	Other works
a In the previous defect models, the oxygen vacancy formation reaction is expressed in terms of (reduction of 2·Mn^4+ to 2·Mn^3+), while the Ereds summarized in Table 2 are for reduction of 2·Mn^3+ to 2·Mn^2+ (as expressed in eqn (1)). Therefore, here we correct the reported reaction energy of the reaction by adding twice of the charge disproportionation reaction energy (i.e. 2·Echg_disp) in the same work to obtain the corresponding Ered.
Oxygen vacancy formation, eqn (1)	E red = 3.43 (eV)	4.51^a (eV)^8
		3.62^a (eV)^10
		2.43^a (eV)^11
		4.11^a (eV)^16
	I red = 2.54 (eV)	—
Cation vacancy formation, eqn (2)	E ox = −4.07 (eV)	−6.54 (eV)^8
		−3.94 (eV)^10
		−2.86 (eV)^11
		−1.73 (eV)^16
	I ox = 16.10 (eV)	—
Charge disproportionation, eqn (3)	E chg_disp = 0.06 (eV)	0.39 (eV)^8
		−1.72 (eV)^11
		0.23 (eV)^16

---

3. The charge disproportionation reaction energy

The LMO charge disproportionation energy calculations are performed at various disproportionated charge (hole and electron) concentrations using different supercell sizes (1 hole/electron in 2 × 2 × 2, 2 × 2 × 4, 2 × 4 × 4, and 4 × 4 × 4 supercells). Following eqn (6), the calculated charge disproportionation energies at different charge concentration are shown in Fig. 5. We extract the charge disproportionation energy of 0.06 eV by taking the linear extrapolated value at the 0 charge concentration, which corresponds to the ideal charge disproportionation reaction energy of an infinitely large LMO supercell. This value is consistent with experimental results, as discussed in Section D.4.

Fig. 5 The calculated ab initio charge disproportionation reaction energies (based on eqn (9)) at different charge concentrations. The extrapolated value at charge concentration = 0 (the black empty circle) is used to represent the charge disproportionation energy in LMO.

4. Ab initio LMO defect energetics vs. the reported defect energetics in literature

Table 2 lists the calculated ab initio LMO defect energetics and the reported values in the literature. The E_rxn parameters correspond to the ideal defect reaction energies of the stoichiometric LMO bulk (the zeroth order defect term), and I_rxn parameters (the coefficients of the first order defect concentration term) correspond to the defect interactions. The ideal oxygen vacancy formation reaction energy (E_red) and cation vacancy formation energy (E_ox) are within the range of the previously reported values. The large negative E_ox suggests formations of cation vacancies is thermodynamically favorable for stoichiometric LMO at the standard condition (T = 298.15 K and P(O₂) = 1 bar). This suggests that at low temperatures and P(O₂) = 1 bar, formation of cation vacancies (oxygen overstoichiometry) will persist until the cation vacancy concentration is high enough that the cation vacancy interactions compensate for the negative cation vacancy formation energy and lead to thermodynamic equilibrium. At high temperatures, the vibrational free energy of the O₂ gas phase will also drive the cation vacancy formation energy in the positive direction (in fact, at 1 atm P(O₂) the cation vacancy formation energy crosses zero and becomes positive at T = 1650 K).

The I_ox of the cation defect vacancy formation reaction (16.10 eV) in Table 2 is consistent with the scale of Mebaneet al.'s strain energy parameters for La_0.9Sr_0.1MnO₃ (22.5 eV), although somewhat less than for La_0.8Sr_0.2MnO₃ (57.9 eV).¹⁵Mebane, et al. reported no value for the LMO system. We hypothesize that the trend of increasing I_ox with increasing Sr content can be understood as additional mean-field energy penalty from Sr-cation vacancy repulsion. From defect modeling perspective, the I_ox parameter in this work is equivalent to the strain field interaction parameter in ref. 15, but our calculated defect interactions also contain other types of defect interactions beyond the strain (e.g. screened electrostatic interactions). A similar form of the defect energy expression (Δε_rxn(c_defect) = E_rxn + I_rxn·c_defect)) can also be found in the rigid band based defect model for electron-rich perovskites, as discussed in ref. 42 and 43. In such case, I_rxn parameters are proportional to be the inverse of the averaged density of states at the Fermi level and accounts for the increase of electron/hole chemical potential with respect to the electron/hole doping. The equations from the rigid band model can be used to estimate density of states at the Fermi level from the I values. However, using these expressions, our normalized I_ox/6 and I_red/2 (normalized as per charge) lead to a inconsistent average density of states at Fermi level (0.787 eV⁻¹ and 0.373 eV⁻¹ per formula unit, respectively). The inconsistency with respect to hole and electron doping suggests that the rigid band type model may not be appropriate to model the LMO/LSM defect chemistry. Overall, the first order dependence of defect formation energy on defect concentration function seems to be an effective defect energy expression for the nonstoichiometric LMO under SOFC conditions as, unlike LSM, this material has no obvious oxygen overstoichiometry flattening/transition in the experimental nonstoichiometry vs.P(O₂) data.

The calculated I_red parameter of the oxygen vacancy formation reaction is only 2.54 eV, which suggests quite effective screening of oxygen vacancy interactions in LMO. This is also why the random defect model (i.e., ideal solution model) generally shows good agreement on the oxygen nonstoichiometry in the oxygen deficient region, and poor fitting is most commonly seen in the oxygen excess region.^11,12

The defect interactions in LSM could arise from strain¹⁵ and/or electronic effects, where the latter may include both electrostatic interactions^44–46 and overall shifts in electron chemical potential.^42,43 To assess the relative contributions of strain and electronic effects we repeated the total energy vs. defect concentration calculations shown in Fig. 4 using ideal (non-distorted) LMO perovskite supercells without relaxation (all atoms in the defected cells were frozen at the undefected LMO atomic coordinates). Because these cells have no relaxation there is no strain interaction between the defects and any interactions are due to electronic contributions. The results of these unrelaxed calculations are provided in Fig. S1 of the Supplemental Information.† The same analysis used on the data from Fig. 4 yields unrelaxed interaction terms (denoted with superscript UR) I^UR_red = 2.66 eV and I^UR_ox = 18.44 eV. These values are quite close to those obtained from the fully relaxed calculations, I_red = 2.54 eV and I_ox = 16.10 eV (see Table 2). The close agreement in values demonstrates that electronic effects dominate LMO defect nonideality. On the other hand, analysis of the first order defect concentration coefficients in Fig. 4 and Fig. S1 demonstrates that relaxations reduce the defect formation energies (E_rxn) by 1.2 eV and 1.5 eV for anion and cation vacancies, respectively. Thus the strain effects do play an important role in the basic defect formation energies. Resolving the relative contributions of electrostatic vs. electron doping effects to the defect interactions is quite difficult due to general challenges in first-principles based supercell calculations.^42–46 Furthermore, the transition in LMO from an insulating system at room temperature to a half metal-like system at SOFC temperatures suggests LMO may behave as a transitional case between a localized (electrostatic interaction) and delocalized (shift in electron chemical potential) electronic system. Resolving the physics of electrons in such system is likely beyond the resolution of the approaches used in this work. Therefore, further studies and improvements on the electronic structure calculation methods are needed to resolve and understand the factors controlling the electronic effect in defect nonideality of LMO.

Finally, we note that the reported charge disproportionation reaction energy varies from −1.74 eV to 0.39 eV (as shown in Table 2) in the previous LMO defect models, which potentially leads to discrepancies in the other two defect energetics. Stevenson et al.⁴⁷ estimated the LMO charge disproportionation reaction free energy based on the measured LMO Seebeck coefficients and the oxygen nonstoichiometry and obtained 0.04 ± 0.03 eV between T = 825 K and 1273 K. The corresponding charge disproportionation enthalpy for the reaction in eqn (3) extracted from the slope of K_D (the equilibrium constant for the charge disproportionation reaction reported in ref. 47) vs. 1/T at T = 825 K and 1273 K is ∼0.05 eV. The close to zero energy for the charge disproportionation reaction is consistent with the closing of the LMO optical gap above T_JT reported by Tobe et al.⁴⁸ Our predicted value of 0.06 eV for the charge disproportion enthalpy is quite consistent with the above experimental measurements. Our calculations, although nominally at zero temperature, make use of cubic structures appropriate to the high-temperature LMO phase. Therefore, our calculated charge disproportionation energy can be compared directly to this measured value.

In the Supplemental Information,† we show the sequential approach, adapted from Poulsen,¹² for solving our LMO defect model with the defect reaction energies expressed by first order defect concentration polynomials. Table 3 summarizes the independent variables to be solved in the LaMnO_3±δ defect model. The resulting stoichiometry predictions will be compared to known experimental results in the next section.

Table 3 Variables to be solved in the LMO defect model

Sublattice	Independent variables	Dependent variables
a Sr doping concentration is specified to be 0 for undoped LMO.
A-site		^a
B-site
O-site

---

E. The predicted results from the defect model

In Fig. 6, the LMO oxygen nonstoichiometry vs.P(O₂) predicted from the above defect model with the calculated ab initio defect energetics are compared with the experimental TG data between T = 873 K and 1273 K reported by Mizusaki et al.¹³ along with the T = 1473 K data from Kamata et al.⁴⁹ These studies provide quite extensive data sets of LMO nonstoichiometry vs.P(O₂) and temperatures. Overall, the results show good agreement with the experimental data, suggesting the calculated ab initio energetics and the estimated vibrational entropies of the defect reactions are capable of describing the defect chemistry of LMO. Not surprisingly, an improved fitting can be further obtained by slightly adjusting the defect energetics and defect interaction parameters. However, such a refitting may not be unique, given the many variables available in the model. In Fig. S2 of the Supplemental Information† we show a specific example of the kind of improved fit which is possible, where we have altered I_red = 2.54 eV from Table 2 to be I_red = 0.8 eV. Furthermore, a slight variation in the oxygen nonstoichiometry vs.P(O₂) is also found in the LMO TG data reported by different groups^8,13,16 and it is not clear that additional fitting to any particularly data set would provide more accurate parameters. In Fig. 7, we show the Brouwer diagram at T = 873, 1073, and 1273 K calculated from the defect model in this work.

Fig. 6 Predictions of LaMnO₃ oxygen nonstoichiometry vs.P(O₂) based on the defect model in this work with the calculated ab initio defect energetics (solid lines). The experimental data are taken from ref. 13 (873 K–1273 K) and ref. 49 (1473 K).

Fig. 7 Calculated Brouwer diagram at T = 873 K, 1073 K, and 1273 K in this work.

The present model helps resolve the large discrepancy (up to 4 eV differences) existing in the LMO defect energies among previous defect models. Furthermore, this work demonstrates the capability of modeling the defect chemistry of LMO with defect energetics directly from ab initio calculations, instead of fitting to experimental LMO oxygen nonstoichiometry.

Since non-ideality in the LMO defect reactions can lead to additional defect concentration dependence in the equilibrium constants, the large error bars of the TG determined equilibrium constants shown in ref. 11 are at least in part a reflection of the non-ideality of defects. In Fig. 8, we show the defect reaction equilibrium constants vs.P(O₂) for the oxygen vacancy formation reaction (K_red) and the cation vacancy formation reaction (K_ox) determined via eqn (S.1) and (S.2) of the Supplemental Information† and the predicted defect concentrations in the Brouwer diagram (Fig. 7). In the case of an ideal solution system, the equilibrium constants should be independent of P(O₂). Such independence can be seen in the intermediate P(O₂) region where cation and oxygen vacancy concentrations are very dilute and the system is close to an ideal solution. However, as the system enters into the non-dilute defect regions, deviations of the equilibrium constants become more significant (specifically, deviations appear for K_ox in the high P(O₂) region and K_red in the low P(O₂) region). These deviations will give rise to errors in the equilibrium constants vs.P(O₂) calculated using the experimental oxygen nonstoichiometry. For example, the equilibrium constant of the cation vacancy formation reaction in LMO can be expressed as:

where δ represents the measured oxygen overstoichiometry in LMO (i.e., LaMnO₃₊δ). Because the TG analysis has limited resolution for small δ, more sampling of δ is generally performed in the high (and low) P(O₂) region, where δ is large and nonideality becomes most significant. Consequently, a large error bar will be found for the derived equilibrium constant based on eqn (10), which then leads to large error bars on the defect reaction enthalpies and entropies. Overall, these results suggest defect energetics in LMO are in fact not constant quantities and that extraction of constant values from the Vant Hoff plots (K_rxnvs. 1/T) based on experimental oxygen nonstoichiometry can lead to significant errors. Care must be taken when interpreting the equilibrium constants obtained from the TG data.

Fig. 8 The defect equilibrium constants (dotted lines for the cation vacancy formation reaction, K_ox, and solid lines for the oxygen vacancy formation reaction, K_red) vs.P(O₂) for LMO at T = 873 K, 1073 K, and 1273 K obtained from the extracted defect concentrations in the Brouwer diagram.

Researchers have also attempted to verify the LSM/LMO defect models by comparing the predicted electrical properties such as Seebeck coefficient and electrical conductivity with experiments^11,12,50 based on specific assumed electron conduction mechanisms. However, these conductivity results seem to be less conclusive.^12,15 We find disagreement in our predicted vs. experimentally measured electronic properties, primarily in the nondilute defect region, which is shown in our estimation of the Seebeck coefficient vs.P(O₂) (via Heikes formula^47,51) and the electronic conductivity vs.P(O₂) fitting (following Poulsen¹²) (data not shown). Although the developed defect model in this work is capable of describing the oxygen nonstoichiometry directly from first principles calculations based on the assumed defect reactions (e.g.eqn (1)–(3)), more work is still needed to further refine the defect model and/or include new physics/energetics in order to explain the electronic properties of LMO/LSM. For example, Kamata et al.⁵² performed electronic conductivity and Seebeck coefficient measurements on LSM suggested that for x ≤ 0.2 both electronic conductivity and Seebeck coefficient can be described in terms of a multi-level hopping conduction model. The multi-level hopping model may suggest the need to include low spin and high spin Mn³⁺ species (with different chemical potentials) in the current LMO/LSM model. We note that the problem in predicting conductivity properties is a limitation shared by the simplified LMO/LSM defect model assumptions. It is therefore a limitation of our present class of defect models and is not a result of using ab initio energetics.

An additional method to verify the defect model is to use predicted oxygen vacancy concentrations in a model of oxygen tracer diffusion and compare to tracer diffusion experiments. The oxygen tracer diffusion can be modeled by the equation:⁵³

where g is the geometry factor (approximated as 1), f₀ is a correlation factor (taken as 0.69 for the oxygen sites in the perovskite lattice⁵⁴), c_vac is the predicted O vacancy concentration of the defect model (), ν₀ is the attempt frequency (taken as 10¹³ s⁻¹), ΔS_m is the migration entropy (taken as 5·k_b), and ΔE_m is the migration energy of oxygen vacancy hopping in LMO. The value of ΔS_m = 5·k_b is set by fitting the value of D* at 1173 K. This value is somewhat large, but may reflect other errors propagating from the rest of the model. The experimental activation barrier of the oxygen tracer diffusion (ΔE_a) in LMO is reported to be ∼2.5 eV by Berenov et al.⁵⁵ in the temperature range of 973–1173 K. A value for ΔE_m is obtained by calculating the oxygen vacancy hopping barrier in a 2 × 2 × 2 supercell with the GGA+U approach (U = 4 eV). We predict ΔE_m is about 0.7 eV, which is consistent with the reported experimental value.⁵⁶ Combining ΔE_m and the predicted oxygen vacancy concentration between 973–1173 K in this work, we obtain an activation barrier of the oxygen tracer diffusion equal to 3.6 eV, which is 1.1 eV higher than the experimental value, as shown in Fig. 9. While such discrepancy seems to be large, the predicted D* using the defect energetic parameters shown in Table 2 are within half an order of magnitude as compared to the experimental D*. Furthermore, it is possible that in the oxygen excess region the cation and anion vacancies have attractive interactions that significantly alter the anion defect formation energies, reducing the effective activation energy for anion diffusion. However, it is also of value to explore if changes in the defect model might be able to explain the observed discrepancy. Therefore, we have performed the exercise of setting experimental D* data of ref. 55 as a constraint and exploring to what extent our ab initio defect energy error bars (∼0.1 eV) could lead to the 1.1 eV difference in activation barrier of D*.

Fig. 9 Predicted of LaMnO₃ oxygen tracer diffusion coefficient (D*) between T = 973 K and 1173 K at P(O₂) = 1 atm based on eqn (11) with the predicted oxygen vacancy concentration in this work and the calculated oxygen vacancy hopping barrier of 0.7 eV. Data A (dashed line) is derived from using a constant charge disproportionation energy E_{chg_disp} = 0.06 eV while data B (solid line) is fitted to the experimental data from ref. 55 with a temperature dependent charge disproportionation energy (E_{chg_disp} = 0.18, 0.10, and 0.02 eV at 973, 1073, and 1173 K).

We observe that a temperature dependent charge disproportionation energy, where E_{chg_disp} = 0.18, 0.10, and 0.02 eV at 973, 1073 and 1173 K, respectively, can be used to bring the model prediction of both the oxygen tracer diffusion activation barrier and the oxygen nonstoichiometry vs.P(O₂) (shown in Fig. 10) in good agreement with the experimental values. While these values are clearly an empirical fit, they show a simple linear decrease with increasing temperature and are quite consistent with experimental understanding. Specifically, the decrease of the charge E_{chg_disp} with increasing temperature is consistent with Goodenough et al.'s vibronic conduction picture of LMO above T_JT, where mobile charge carriers are strongly coupled to the bond-length fluctuations, and hence the localized electrons in LMO are approaching the transition from localized to itinerant electronic behavior with increasing temperature.⁵⁷ In addition, all the fitted E_{chg_disp}s between 973–1173 K are smaller than the reported experimental optical gap of 0.2 eV measured at 800 K.⁴⁸ Therefore, we take as a simple semi-empirical defect model that the charge disproportionation energy is fairly constant at 0.2 eV up to 948 K, at which point it reduces linearly to 0 (itinerant behavior) at 1198 K. This model for charge disproportionation is consistent with experimental data on the gap and known oxygen diffusion data. Furthermore, using this model, the predicted oxygen nonstoichiometry vs.P(O₂) between 873 K to 1473 K is still in excellent agreement with experiments, as shown in Fig. 10. We can see from these results that oxygen vacancy concentration at high P(O₂) is very sensitive to the charge disproportionation energy. Such a sensitivity is perhaps beyond the capability of present ab initio models, but future improvements will hopefully remove the need for such additional fitting. Finally, we note that even using the modified E_{chg_disp}(T), the discrepancy between the predicted electronic conductivity vs. experiments still persists, since the hole carriers are not sensitive to the small change in E_{chg_disp}(T). Hence, the defect model still needs to be refined to include new physics/energetics in order to explain the electronic properties of LMO. Nonetheless, to our knowledge, this work presents the first LMO defect model that agrees with both the defect thermodynamic and the oxygen tracer diffusivity measurements.

Fig. 10 Predicted of LaMnO₃ oxygen nonstoichiometry vs.P(O₂) between T = 873 K and 1473 K based on the defect model in this work along with modified temperature dependent charge disproportionation energies (see text) and the calculated ab initiooxygen vacancy and cation vacancy formation defect energetics in Table 2 (solid lines). The experimental data are taken from ref. 13 and 49.

F. Conclusions

An ab initio and empirical defect model for LMO under SOFC conditions has been developed. The model explicitly calculates the reaction enthalpies and entropies with ab initio methods, although empirical values for O₂ gas are used and some fitting to experimental data is needed to predict the oxygen diffusivity. A useful outcome of analysing this model is the realization that the estimated vibrational entropy of the defect reactions in LMO involving O²⁻(solid) ↔ ½·O₂(gas) + 2·e⁻ is about ∼80 J/(mol K) at T = 800–1300 K and relatively insensitive to the specific details of the solid phase vibrational model. This provides a useful constraint on parameter values for experimental fits of defect models. The calculated ab initio defect energetics at various defect concentrations are expressed with first order functions of defect concentration to describe nonideality of defects in terms of solution thermodynamics. Defect interactions have been challenging to include in LMO defect models fit to experiment as the interaction parameters are difficult to constrain from thermodynamic data. The defect interactions are predicted to be quite significant, particularly for the cation vacancies. The influence of these interactions at least partially explains the range of defect reaction energies that have been found in fitting models that do not explicitly allow for defect interactions. The predicted oxygen nonstoichiometry vs.P(O₂) from the model shows good agreement with the experimental data over a wide range of temperatures (T = 873 K–1473 K) and P(O₂), suggesting that both the ab initio defect energetics and the vibrational entropies of defect reactions we have obtained are capable of describing the defect chemistry of LMO. Previous LMO/LSM models have defect energies that vary by many electron volts. This work provides a set of ab initio based defect reaction energies that can be used for accurate modeling of the LMO system, greatly reducing the uncertainty in the key defect reaction energies for this material. The predicted oxygen tracer diffusion kinetics based on the defect concentration derived from the defect model with a constant charge disproportionation energy of 0.06 eV shows 1.1 eV overestimation on the oxygen tracer diffusion activation barrier as compared to the experimental value. However, by incorporating an empirical temperature dependent charge disproportionation energy decreasing linearly from 0.2 eV to 0 eV from 948 K to 1198 K, the predicted defect concentration from the defect model is capable of describing both the experimental oxygen tracer diffusion kinetics and defect thermodynamics in LMO. Such temperature dependence for the charge disproportionation reaction is consistent with the vibronic conduction picture of LMO above T_JT as proposed by Goodenough et al.,⁵⁷ and the magnitude is within the reported experimental optical gap (0.2 eV) above T_JT. Overall, this work depicts a framework for modeling defect thermodynamics initiated from ab initio calculations with experimental data inputs as constraints to resolve and understand defect chemistry of complex oxide systems. Similar approaches can be taken to model defect chemistry in a wide range of complex oxides.

Acknowledgements

We gratefully acknowledge financial support from the National Science Foundation MRSEC program under grant number DMR-0520527 and the U.S. Department of Energy (U.S. DOE), Office of Basic Energy Sciences, Division of Materials Sciences and Engineering (under Award No. DE-SC0001284). We also gratefully acknowledge computing support from National Science Foundation (NSF) National Center for Supercomputing Applications (NCSA) under award number DMR060007 and NERSC allocation of the Center for Nanophase Materials Sciences (CNMS) at Oak Ridge National Laboratory under grant number CNMS2008-204.

References

1. S. C. Singhal, High-temperature Solid Oxide Fuel Cells: Fundamentals, Design and Applications, Elsevier Science, 2004.
2. S. B. Adler, Chem. Rev., 2004, 104, 4791–4843.
3. S. McIntosh, J. F. Vente, W. G. Haije, D. H. A. Blank and H. J. M. Bouwmeester, Solid State Ionics, 2006, 177, 1737–1742.
4. J. A. Kilner, R. A. De Souza and I. C. Fullarton, Solid State Ionics, 1996, 86–88, 703–709.
5. R. A. De Souza and J. A. Kilner, Solid State Ionics, 1999, 126, 153–161.
6. R. Merkle, J. Maier and H. J. M. Bouwmeester, Angew. Chem., Int. Ed., 2004, 43, 5069–5073.
7. S. B. Adler, J. A. Lane and B. C. H. Steele, J. Electrochem. Soc., 1996, 143, 3554–3564.
8. J. H. Kuo, H. U. Anderson and D. M. Sparlin, J. Solid State Chem., 1989, 83, 52–60.
9. J. A. M. Vanroosmalen and E. H. P. Cordfunke, J. Solid State Chem., 1991, 93, 212–219.
10. J. A. M. Vanroosmalen and E. H. P. Cordfunke, J. Solid State Chem., 1994, 110, 113–117.
11. J. Nowotny and M. Rekas, J. Am. Ceram. Soc., 1998, 81, 67–80.
12. F. W. Poulsen, Solid State Ionics, 2000, 129, 145–162.
13. J. Mizusaki, N. Mori, H. Takai, Y. Yonemura, H. Minamiue, H. Tagawa, M. Dokiya, H. Inaba, K. Naraya, T. Sasamoto and T. Hashimoto, Solid State Ionics, 2000, 129, 163–177.
14. K. Nakamura, J. Solid State Chem., 2003, 173, 299–308.
15. D. S. Mebane, Y. J. Liu and M. L. Liu, Solid State Ionics, 2008, 178, 1950–1957.
16. A. Y. Zuev and D. S. Tsvetkov, Solid State Ionics, 2010, 181, 557–563.
17. L. Malavasi, C. Ritter, M. Cristina Mozzati, C. Tealdi, M. Saiful Islam, C. Bruno Azzoni and G. Flor, J. Solid State Chem., 2005, 178, 2042–2049.
18. J. A. M. Vanroosmalen, E. H. P. Cordfunke, R. B. Helmholdt and H. W. Zandbergen, J. Solid State Chem., 1994, 110, 100–105.
19. J. F. Mitchell, D. N. Argyriou, C. D. Potter, D. G. Hinks, J. D. Jorgensen and S. D. Bader, Phys. Rev. B: Condens. Matter, 1996, 54, 6172–6183.
20. H. L. Lukas, S. G. Fries and B. Sundman, Computational Thermodynamics: The Calphad Method, Cambridge University Press, Cambridge, UK, 2007.
21. Y. A. Chang and W. A. Oates, Materials Thermodynamics, John Wiley & Sons Inc., Hoboken, New Jersey, USA, 2009.
22. NIST, in NIST Standard Reference Database No. 69, ed. P. J. Linstrom and W. G. Mallard, National Institute of Standards and Technology, Gaithersburg, MD, 2003, http://webbook.nist.gov/chemistry/.
23. K. Reuter and M. Scheffler, Phys. Rev. B: Condens. Matter, 2001, 65, 035406.
24. E. A. Kotomin, Y. A. Mastrikov, E. Heifets and J. Maier, Phys. Chem. Chem. Phys., 2008, 10, 4644–4649.
25. Y.-L. Lee, J. Kleis, J. Rossmeisl and D. Morgan, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 80, 224101.
26. Y. A. Mastrikov, R. Merkle, E. Heifets, E. A. Kotomin and J. Maier, J. Phys. Chem. C, 2010, 114, 3017–3027.
27. J. Callaway, Quantum Theory of the Solid State, Academic Press, 1974.
28. M. N. Iliev, M. V. Abrashev, H. G. Lee, V. N. Popov, Y. Y. Sun, C. Thomsen, R. L. Meng and C. W. Chu, Phys. Rev. B: Condens. Matter Mater. Phys., 1998, 57, 2872.
29. F. Reif, Fundamentals of Statistical and Thermal Physics (McGraw-Hill Series in Fundamentals of Physics), McGraw-Hill Science/Engineering/Math, 1965.
30. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter, 1993, 47, 558.
31. G. Kresse and J. Furthmuller, Phys. Rev. B: Condens. Matter, 1996, 54, 11169–11186.
32. P. E. Blochl, Phys. Rev. B: Condens. Matter, 1994, 50, 17953–17979.
33. G. Kresse and D. Joubert, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 1758–1775.
34. J. P. Perdew and Y. Wang, Phys. Rev. B: Condens. Matter, 1992, 45, 13244.
35. V. I. Anisimov, J. Zaanen and O. K. Andersen, Phys. Rev. B: Condens. Matter, 1991, 44, 943–954.
36. V. I. Anisimov, F. Aryasetiawan and A. I. Lichtenstein, J. Phys.: Condens. Matter, 1997, 9, 767–808.
37. L. Wang, T. Maxisch and G. Ceder, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 73, 195107.
38. S. Fritsch and A. Navrotsky, J. Am. Ceram. Soc., 1996, 79, 1761–1768.
39. E. Heifets, J. Ho and B. Merinov, Phys. Rev. B: Condens. Matter Mater. Phys., 2007, 75, 155431.
40. S. B. Adler, J. Am. Ceram. Soc., 2001, 84, 2117–2119.
41. X. Y. Chen, J. S. Yu and S. B. Adler, Chem. Mater., 2005, 17, 4537–4546.
42. M. H. R. Lankhorst, H. J. M. Bouwmeester and H. Verweij, Phys. Rev. Lett., 1996, 77, 2989–2992.
43. M. H. R. Lankhorst, H. J. M. Bouwmeester and H. Verweij, J. Am. Ceram. Soc., 1997, 80, 2175–2198.
44. G. Makov and M. C. Payne, Phys. Rev. B: Condens. Matter, 1995, 51, 4014–4022.
45. C. W. M. Castleton, et al. , Modell. Simul. Mater. Sci. Eng., 2009, 17, 084003.
46. R. M. Nieminen, Modell. Simul. Mater. Sci. Eng., 2009, 17, 084001.
47. J. W. Stevenson, M. M. Nasrallah, H. U. Anderson and D. M. Sparlin, J. Solid State Chem., 1993, 102, 175–184.
48. K. Tobe, T. Kimura, Y. Okimoto and Y. Tokura, Phys. Rev. B: Condens. Matter, 2001, 64, 184421.
49. K. Kamata, T. Nakajima, T. Hayashi and T. Nakamura, Mater. Res. Bull., 1978, 13, 49–54.
50. J. H. Kuo, H. U. Anderson and D. M. Sparlin, J. Solid State Chem., 1990, 87, 55–63.
51. R. R. Heikes, in Thermoelectricity: Science and Engineering, ed. Robert R. Heikes and Roland W. Ure, Jr., with the collaboration of Stephen J. Angello and others, Interscience Publishers, New York, 1961.
52. H. Kamata, Y. Yonemura, J. Mizusaki, H. Tagawa, K. Naraya and T. Sasamoto, J. Phys. Chem. Solids, 1995, 56, 943–950.
53. H. Mehrer, Diffusion in Solids, Springer-Verlag, New York, 2007.
54. T. Ishigaki, S. Yamauchi, K. Kishio, J. Mizusaki and K. Fueki, J. Solid State Chem., 1988, 73, 179–187.
55. A. V. Berenov, J. L. MacManus-Driscoll and J. A. Kilner, Solid State Ionics, 1999, 122, 41–49.
56. A. Belzner, T. M. Gür and R. A. Huggins, Solid State Ionics, 1992, 57, 327–337.
57. J. B. Goodenough, J. S. Zhou, F. Rivadulla and E. Winkler, J. Solid State Chem., 2003, 175, 116–123.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c1cp22380a

---