First-principles study of the morphology and surface structure of LaCoO3 and La0.5Sr0.5Fe0.5Co0.5O3 perovskites as air electrodes for solid oxide fuel cells

ABSTRACT The slow kinetics of the oxygen exchange reaction at the surface of air electrodes for solid oxide fuel cells must be improved to reduce their operating temperature. In this study, we use first-principles computation, based on the density functional theory, to investigate the surface morphology and oxygen vacancy formation energies of perovskite-type LaCoO3 and La0.5Sr0.5Co0.5Fe0.5O3 (LSCF5555) compounds. The cationic arrangement is optimized energetically for LSCF5555 using a genetic algorithm approach. We confirm that Sr and Co condensation occurs near the LSCF5555 surface. Significant lowering of the oxygen vacancy formation energies at the perovskite surface is indicated for both the compounds. In addition, several surface O sites in LSCF5555 exhibit lower oxygen vacancy formation energies than those in LaCoO3, which may improve the kinetics of the surface charge-transfer reaction in the LSCF-related compounds. Graphical Abstract


Introduction
Perovskite oxide -LaCoO 3-δ -and its derivatives are being investigated as cathode electrode materials for solid oxide fuel cells (SOFCs), owing to their rapid mixed oxide ion and electron conductivities [1][2][3][4]. Using materials with high mixed conductivity enables a reduction in the operating temperature of SOFCs to a moderate range (i.e. 500-700 °C). This approach increases the long-term stability of the device by minimizing the rate of degradation reactions that limit the fuel cell performance. The cathode reaction involves the adsorption of O 2 molecules, followed by an electron exchange reaction to form O 2at the particle surfaces and finally the migration of O 2and electrons in the bulk of the particles. In previous studies, Kawada et al. [5,6]. and Adler et al. [7]. reported that the surface process is one of the rate-determining steps for ion exchange, although the rate-determining step generally depends on the particle size and morphology. Therefore, material optimization of LaCoO 3based compounds, in terms of the surface structure, is necessary for their application in fuel cells.
To date, researchers have focused on studying the defect chemistry and oxygen transport in LaCoO 3 and related materials. For example, the defect chemistry and oxide ion migration in the bulk of these compounds have been studied using transient thermogravimetry (TG) [8,9], coulometric titrations [10,11], and electrical conductivity measurements [12][13][14]. Furthermore, the mechanisms responsible for the above-mentioned properties have been analyzed using quantum beam techniques, such as diffraction and spectroscopy [15,16], and computational techniques [17][18][19]. Similarly, the kinetics of oxygen exchange at the surface have been experimentally studied [20,21], whereas the atomic or electronic mechanisms of the oxygen exchange reaction, at the surface of LaCoO 3 -derived compounds, remain unclear. Recent studies involving low-energy ion scattering techniques have reported that the surface chemistries, i.e. compositions of Sr-doped LaCoO 3 and its derivatives play an important role in the surface reaction kinetics [22][23][24]. Cao et al. demonstrated the critical role of surface defect formation in the kinetics of the surface O exchange reaction through density functional theory (DFT) calculations [20]. Thus, understanding the surface structures of LaCoO 3based compounds is necessary to effectively design their surface compositions.
Hence, in this study, we investigated the surface structures and morphologies of LaCoO 3 and its derivative, La 0.5 Sr 0.5 Co 0.5 Fe 0.5 O 3 (LSCF5555), using firstprinciples DFT calculations. In particular, determining the energetically optimized cation arrangement in the bulk of LSCF5555 is crucial for understanding the aforementioned surface reaction kinetics. However, optimizing the cation arrangement in LSCF5555 is technically difficult, owing to its huge configurational freedoms and limited computational resources. Moreover, despite the importance of surface chemistry in functional materials, only a few attempts have been made to optimize the surface and/or interface structures for nonstoichiometric crystalline solids, such as solid solution LSCF5555. Although DFT-derived Monte Carlo simulations involving the cluster expansion (CE) approach are often used for optimizing cation arrangements [25], the formalism of the CE approach requires that ionic arrangements in crystals possess relatively high symmetry. In this respect, the CE approach is inappropriate as the termination of the periodicity in the surface models of crystalline solids leads to symmetry breakages. Ab initio structure prediction techniques, such as USPEX [26,27] and CALYPSO [28][29][30], are alternatives for developing structural models, including interfaces. These methods are useful, even if an a priori structural framework cannot be obtained for a given composition. However, the ab initio structure prediction methods become computationally demanding when the number of atoms increases and the composition becomes complicated. In particular, nonstoichiometry, such as by partial substitution of atoms, leads to combinatorial explosion, which makes ab initio structure optimization highly expensive in terms of computational resources and time. Instead, a genetic algorithm (GA) based on the specific crystal structure (perovskite) model is used in this study to efficiently determine a surface cation arrangement within a reasonable computational cost. Although the technique requires an a priori structural model using USPEX and CALYPSO, the computational cost reduces significantly by eliminating high structure flexibility. In addition, the crystal and electronic structures are discussed to improve our understanding of the reaction kinetics at the surface of LaCoO 3 and its derivatives.

First-principles density functional theory calculations
All first-principles DFT calculations, for LaCoO 3 , LSCF5555, and their related structures, were performed using the generalized gradient approximation for solids (GGA-PBEsol) [31,32] and GGA + U [33]. We employed exchange-correlation spin-polarized functionals with the plane-wave basis set and projector-augmented wave method [34] and executed the calculations using the Vienna Ab-initio Simulation Package (VASP) [35,36]. The effective U values for the Co and Fe d electrons were set as 3.0 and 4.0 eV, respectively, based on the previous reports [37,38]. The ferromagnetic configurations for the transition metal ions were set as the initial inputs in this study to avoid using large number of spin up/down configurations for the surface model calculations and cation arrangement optimization. A kinetic cutoff energy of 500 eV and k-point meshes were set; the product of meshes and atoms was approximately 1000, which was determined through a convergence test (<5 meV/ABO 3 ). Integration over the first Brillouin zone was performed with Gaussian smearing (σ = 0.05 eV). Structural relaxation, including lattice parameters and fractional coordinates, was performed for the 80-atom rhombohedral LaCoO 3 bulk structure (a 2 × 2 × 2 supercell of the unit cell) according to the previous reports [39,40].
The surface structures of LaCoO 3 were simulated using the slab technique, in which sets of infinite layers, separated by vacuum layers (10 Å in thickness), were repeated periodically along the surface normal [41][42][43]. The slabs were modeled to exhibit inversion or mirror-type symmetry with respect to the center of the slab. Five low-index facets of the {001}, {110}, {111}, {1-10}, and {1-11} surfaces, with various termination structures, were computationally modeled. The stoichiometric and nonstoichiometric O slabs were considered in this study, whereas the cation compositions remained unchanged. Charge-neutral and stoichiometric slabs were constructed by modifying the number of surface O atoms, under the abovementioned symmetry constraints. The surface energy (surface tension), ϒ, for the O nonstoichiometric slab model, was determined using the following equation: where E slab and E bulk refer to the total energies of the slab and the bulk of perovskite compounds, respectively; n and S are the number of moles of the perovskite phase and the surface area, respectively; μ i is the chemical potential of species i; and N i is the number of moles of ions deviating from the stoichiometric composition. In this study, we considered the chemical potential as a temperature-dependent function according to the following equation: where μ � i , R, T, and a i correspond to the chemical potential of species i under standard conditions, gas constant, temperature, and activity of species i, respectively. As we consider O to be nonstoichiometric, activity a i can refer to the O partial pressure, p(O 2 ). Further details are described in [44]. In this study, the other temperature-dependent contributions, such as the vibrational enthalpy and entropy or configurational entropy, were ignored. The thermodynamic equilibrium morphology of a crystal was determined using the Wulff construction procedure [45].

Genetic algorithm optimization for surface structure
Building a suitable configurational model for the bulk and surface of LSCF5555 is a technical challenge. Monte Carlo techniques are often utilized to counter the above-mentioned disadvantages [46]; however, the parameterization of interactions is always challenging, particularly for a slab model. In this study, the configuration was determined via a GA approach by optimizing the La/Sr (at Perovskite-A site) and Co/Fe (at Perovskite-B site) arrangements, where the GA did not require the interaction parameters (all the calculations were performed using DFT approaches), which is in-line with our previous reports [47][48][49][50][51][52]. The exact computational compositions, used for the GA approach, were La 8 Sr 8 Co 8 Fe 8 O 48 and La 10 Sr 10 Co 10 Fe 10 O 60 for the bulk, and the {110} and {1-10} surface models of LSCF5555, respectively. The total number of combinations of cation arrangements exceeded 10 8 or 10 10 for the bulk and surface models in this study. The cation arrangement at Perovskite-A sites was considered as a binary string, consisting of two labels, namely 0 (La) and 1 (Sr), and each cationic site was assigned to a specific index of the string. Likewise, a second binary string, consisting of 0 (Co) and 1 (Fe), was considered for Perovskite-B sites. To preserve the composition, we used a permutation encoding technique for the above binary string representations. We first prepared 20 configurations (binary strings), where La/Sr and Co/Fe were arranged randomly at the cation sites (1st generation), and computed the total electron energies using DFT. Among the 20 configurations, 12 low-energy configurations (60%) were selected as survivors, and their structural features were passed to the next generation using the following four options: (1) three best structures succeeded as they were, (2) eight new structures were created using the two-point crossover technique, (3) eight were created using the halfuniform crossover technique, and (4) four were created using mutation techniques. The total energies for the thus obtained 20 new configurations, using options (2)-(4), were calculated. Twelve survivors were again evaluated, including three best structures. This scheme (generation) was repeated and the lowest-energy configuration was determined heuristically. We obtained the DFT-total energies during the GA procedure by reducing the energy cutoff and number of k points in the reciprocal cell to 360 eV and unity (Γ point sampling), respectively, to speed up each calculation because the total energies for more than 1000 configurations per composition had to be calculated. After the GA cycles, the obtained five lowest-energy configurations were recalculated with sufficient energy cutoff and k-point meshes, and the lowest-energy configuration was chosen. In-house software and VASP were used for GA and DFT calculations, respectively.

Results and discussion
Three types of spin states, namely, high, low, and intermediate, were calculated for Co ions in the bulk rhombohedral LaCoO 3 phase [53]. We could confirm that all the results converged to the input ferromagnetic spin configurations. We found that the electronic structure is metallic for the intermediate-spin state, whereas it is insulating for the low-or high-spin states, agreeing with the findings of a previous study [32]. The intermediate-spin state at 0.23 or 0.42 eV/LaCoO 3 is more stable than the low-or high-spin states in ferromagnetic spin configurations, respectively. Therefore, the intermediate-and low-spin states (1 or 0 μ B as the local spin moment for a Co ion) are considered, hereafter, as the initial spin states of Co ions. Table 1 lists the calculated surface energies for the stoichiometric LaCoO 3 . Lower surface energies are obtained for the intermediate-spin configuration. The results agree with those of a previous study [54], as indicated by the lowest-surface-energy facet {110}, which, in a rhombohedral symmetry, corresponds to the {100} facet in cubic perovskite symmetry. Moreover, some of the calculated surface energies are slightly lower than those given in a previous study [55], which may be due to differences in the calculation settings or input structure model. In principle, DFT calculations seek the most stable electronic/crystal structures. As low-spin states are energetically metastable, electronic structures may partly change into low-energy spin states, depending on the input structure and calculation conditions. To avoid the risk of low reproducibility, we used the surface energies corresponding to the intermediate-spin states, listed in Table 1, for further calculations, unless mentioned otherwise. According to Wulff's theorem [45], the crystal morphology at equilibrium is predicted to be minimal for a given volume, from the surface energies. Figure 1 shows the optimized crystal morphology Nonstoichiometric surface structures with oxygen, such as LaCoO 3±δ , were also modeled to consider the effect of changes in the oxygen chemical potential, such as those introduced by changes in the partial pressure and temperature of oxygen that occur during the operation of an SOFC. In this study, only three facets, that is, {110} (LaO terminated), {1-10} (O 2 terminated), and {100} (LaO 3 terminated), were considered as the major surfaces, as indicated by the Wulff shape of the stoichiometric LaCoO 3 (Figure 1). The surface energies for the nonstoichiometric oxygen model were also calculated as a function of the chemical potential of oxygen, that is, O 2 partial pressure. The molar ratio of the Perovskite-A and B cations was fixed at 1:1 in this study. The total energy of the O 2 molecule was corrected by the experimental s-or p-block metal oxidation reaction energies at an oxygen partial pressure, p(O 2 ), of 0.2 atm, (1 atm total pressure) and temperature of 298 K [56]. The temperature dependence of enthalpy and entropy for the O 2 molecule was taken from previous experimental datasets [57]. Figure 2 (a) shows the temperature dependence of the surface energies of the three facets, where positive and negative changes in the surface energy versus temperature indicate regions with oxygen excess and oxygen deficiency, respectively, whereas no changes indicate that the stoichiometry of the surface composition is maintained. Over the temperature range of 400-2500 K, the LaO-terminated {110} facet is the most stable one. Figure 2 (b) and 2 (c) show the ratio of the surface area for each facet and the morphology as a function of temperature, respectively, based on the calculated surface energies and Wulff's theorem. The results show that the change in morphology is marginal. It should be noted that the present computations only considered the chemical potential of oxygen as the temperature-dependent term. Hence, other factors, such as spin-state-derived phase transitions [58], should be considered in the future studies. Figure 3 (a) shows an example of a GA-driven structure optimization for La/Sr and Co/Fe arrangements in bulk LSCF5555. The total energy distribution of GA-derived structures is reduced as the number of GA generations increased. The best structure (lowest energy) is shown in Figure 3 (b), where no clear cation   ordering is observed, which is in good agreement with experimental results [59]. In addition, the GAoptimized cationic arrangements at the surface of LSCF5555 are shown in Figure 4. A relatively small cell for surface calculations (<100 atoms) was used due to the limited computational resources (more than 3500 DFT calculations were performed in the present GA). In both cases, Co and Sr ions tended to be concentrated at the surface of LSCF5555. Alkaline earth termination is confirmed in the perovskite-type materials according to the experimental observations [22]. Sr ion condensation at the surface of (La 1-x Sr x ) 2 CoO 4 with a Ruddlesden-Popper structure has also been reported in recent experimental results [24]. We also observed Co condensation near the surface of LSCF5555 via transmission electron microscopy analysis, as described in Supplementary Figure S1 and Supplementary Section S1.
To clarify the role of the Sr and Co concentrations at the surface of LSCF5555 in the redox reaction properties, we calculated the oxygen vacancy formation energies at the surface, which are related to the diffusion of oxygen and the formation of oxygen adsorption/desorption sites. The following equation was used to evaluate the oxygen vacancy formation energies, E v : where E def and E st correspond to the total electron energies for defective and stoichiometric compounds, respectively. A single oxygen vacancy was introduced at various oxygen sites in the stoichiometric slab model for E def calculations. Here, μ O 2 is the chemical potential of oxygen, for which we used the total energy of the O 2 molecule in a vacuum. Figures 5 and 6 show the calculated oxygen vacancy formation energies for the {110} and {1-10} facets of LaCoO 3 or LSCF5555 perovskites, respectively, as a function of the z-fractional coordinates of the oxygen atom in the slab model. The oxygen vacancy formation energies are scattered even for the same z-fractional coordinates of the LaCoO 3 compounds owing to the symmetry breakage due to the termination of the periodicity at the crystal surfaces and minor numerical errors. The average and standard deviations of the oxygen vacancy formation energies as a function of the z-fractional coordinates are displayed in Supplementary Figure S2. Around the center of the {110} facet slab model ( Figure 5), the calculated oxygen vacancy formation energies for LSCF5555 were slightly smaller than or comparable to those of LaCoO 3 , where the energies were in the range of 3.0-3.5 eV. Moreover, the oxygen vacancy formation energies were scattered for LaCoO 3 , and decreased to approximately 2.0 eV for LSCF5555 close to the surface. Therefore, a higher concentration of O vacancies was observed at the surface of LSCF5555 than that of the bulk of LSCF5555 or LaCoO 3 . A similar trend was observed for the {11-0} facet slab model, in which a lower oxygen vacancy formation energy was clearly visible at the surface than that of the bulk of either LaCoO 3 or LSCF5555. The vacancy formation energies at the surface for LSCF5555 were slightly smaller than those for LaCoO 3 . Hence, we inferred that a higher oxygen vacancy concentration at the surface of LSCF5555 may have been a source of the observed higher catalytic activity and/or electrochemical performance [22]. Parts (b) and (c) of Figures 5 and 6 show the Bader charge distribution of Co and oxide ions, which can be used to understand the distribution of oxygen vacancy formation energies in the slab model of LaCoO 3 and/or LSCF5555 compounds. The Bader charge is determined by the partitioning of the electronic density according to its zero-flux surfaces [60][61][62]. We removed the plots of the Bader charges for the ions located at the surface, as a lack of surrounding ions results in an overestimation of the number of electrons. The distributions of the Bader charges for the Co and oxide ions were greater close to the surface than to the bulk, that is, the center of the slab model, for both compounds. In particular, the oxide ions with higher Bader charges exhibited low oxygen vacancy formation energies, whereas no significant relationship was observed for the Bader charges of the surrounding Co ions and the oxygen vacancy formation energies. Hence, Co condensation at the surface is not the primary reason for the observed increase in the Bader charge and decrease in the vacancy formation energies of oxide ions at the surface. We can infer that lowervalency Sr 2+ condensation enhances the increase in the Bader charge of oxide ions at the surface due to charge compensation.

Conclusions
The surface structures and oxygen vacancy formation energies of perovskite-type compounds, LaCoO 3 and LSCF5555, as candidate cathode materials for SOFC, were investigated using DFT-GGA+U computations. The surface energy calculations for the low-index facets of LaCoO 3±x indicated that the {110} and {1-10} facets were the dominant surfaces for the Wulff-theoremderived crystal shapes. In addition, the GA optimization of the arrangement of ions revealed that Sr and Co ions are concentrated at both facet models of the LSCF5555 compounds. The lowering of the oxygen vacancy formation energy was observed at the surface of both compounds compared with that in the bulk. Specifically, some oxide ion sites at the surface of LSCF5555 exhibited lower vacancy formation energies than that of LaCoO 3 . This explains the better SOFC performance when using LSCF than LaCoO 3 .

Disclosure statement
There are no conflicts of interest to declare.

Data availability statement
The data that support the findings of this study are available from the corresponding author upon reasonable request.

Funding
This work was partially supported by the Grant-in-Aid for Scientific Research under Grants 18K19129, 18H01707, 19H05815, and 20H02436 from the Ministry of Education Culture, Sports, Science and Technology, Japan (MEXT); the Elements Strategy Initiative to Form Core Research