A computational study of the electronic properties, ionic conduction, and thermal expansion of Sm_1−xA_xCoO₃ and Sm_1−xA_xCoO_3−x/2 (A = Ba²⁺, Ca²⁺, Sr²⁺, and x = 0.25, 0.5) as intermediate temperature SOFC cathodes

Abstract

The substitutional doping of Ca²⁺, Sr²⁺, and Ba²⁺ on the Sm-site in the cubic perovskite SmCoO₃ is reported to improve both electronic and ionic conductivities for applications as solid oxide fuel cell (SOFC) cathodes. Hence, in this study we have used density functional theory (DFT) calculations to investigate dopant configurations at two different dopant concentrations: 25 and 50%. To preserve the electroneutrality of the system, we have studied two different charge compensation mechanisms: the creation of oxygen vacancies, and electronic holes. After examining the electronic structure, charge density difference, and oxygen vacancy formation energies, we concluded that oxygen vacancy charge compensation is the preferred mechanism to maintain the electroneutrality of the system. Furthermore, we found that the improvement of the electronic conduction is not a direct consequence of the appearance of electron holes, but a result of the distortion of the material, more specifically, the distortion of the Co–O bonds. Finally, molecular dynamics were employed to model ionic conduction and thermal expansion coefficients. It was found that all dopants at both concentrations showed high ionic conduction comparable to experimental results.

---

1. Introduction

Solid oxide fuel cells (SOFC) are important potential replacements of traditional power sources, and offer a clean and efficient chemical-to-electrical energy conversion process.^1–3 However, to make these devices more widely available and applicable, their operating temperature and costs have to be lowered. One of the most pressing issues is the need to develop new cathode materials that are capable of operating at intermediate temperatures (500–900 °C), with high electrical and catalytic efficiencies at these lower temperatures, i.e. intermediate temperature SOFC (IT-SOFC). The factor influencing the efficiency of IT-SOFC cathodes is the catalytic efficiency towards the oxygen reduction reaction, which depends on the surface oxygen reduction and oxygen bulk diffusion. These reactions are oxygen vacancy assisted, whereby high oxygen vacancy concentrations are favorable. The main task of the IT-SOFC cathode is then to provide an efficient pathway for the reduced oxygen through the cathode bulk to the electrolyte, as well as being able to transport electrons.¹ Therefore, IT-SOFC cathodes are required to have high ambipolar conductivity in the form of mixed electronic and ionic conductivity (MIEC).⁴ Furthermore, it is critical that the cathode is chemically compatible with the electrolyte, and that all components have similar thermal expansion coefficients, so as to have physical compatibility.

Unlike traditional cathode materials such as La_1−xSr_xMnO_3−δ (LSM), which loses efficiency at these low temperatures due to polarization resistance, cobalt-based perovskites have shown better performance and have been suggested as suitable materials for IT-SOFC,⁵ as their reported properties include better mixed ionic and electronic conduction, low polarization resistance, and lower overpotentials.^6,7 Their higher ionic conductivity is a consequence of the lower activation energies required for the oxygen to migrate via oxygen vacancies, which are induced when substituting Sm³⁺ by a divalent cation (A).^6–13 Upon doping, the overall charge of the system can be compensated via two different charge compensation schemes. On the one hand, oxygen vacancies (V_O) can be generated, as described in eqn (1), according to the Kröger–Vink notation:¹⁴

On the other hand, the charge can be compensated with the creation of electron holes, as expressed in eqn (2).^13,15–17

While the oxygen vacancy compensation mechanism is widely used to explain the formation of oxygen vacancies, and thus the increase in the ionic conductivity, the hole compensation mechanism has been suggested to explain the increase in the electronic conductivity.^16,17 Although there is no strong evidence, it is believed that both mechanisms coexist in the material. The increase in the electronic conductivity seems to be related with the spin transition observed in cobalt, from low spin to either high (HS, t⁴_2ge²_g), intermediate (IS, t⁵_2ge¹_g), or mixed spin (LS/IS, IS/HS).⁸ This leads to magnetic domains, metallic electronic structure, and the cobalt ions coexisting in different spin states.⁸ Another improvement that doped SmCoO₃ cathodes present with respect to LSMO, is that although they do not couple well with yttria-stabilized zirconia (YSZ) as the electrolyte material, they show better compatibility with ceria-based electrolytes such as gadolinium-doped ceria (GDC).^6,10

Although experimental studies conducted on Sm_1−xA_xCoO_3−x/2 indicate that they could be promising materials for IT-SOFC cathodes, questions remain why samarium cobaltates show better performance than doped lanthanum manganite.² Hence, in this paper, we investigate three common SmCoO₃ dopants; Ca²⁺, Sr²⁺, and Ba²⁺, at two different dopant concentrations (x); x = 0.25, and 0.5. Using a combination of density functional theory and molecular dynamics simulations, we have investigated two different charge compensation schemes, the electronic and magnetic structure, oxygen diffusion, and the thermal expansion coefficients. On the basis of these calculations, we then compare the different Sm_1−xA_xCoO_3−x/2 systems to gain an understanding of how the dopant species and its concentration influences the oxygen conduction and the materials’ properties. The results have been compared to experiment, and especially the compatibility with traditional electrolytes has been kept in mind when evaluating which dopant and concentration are most suitable for IT-SOFC cathode applications.

2. Computational details

In this study, we have employed two different methodologies to model the Sm_1−xA_xCoO_3−x/2 systems; density functional theory (DFT) calculations and interatomic potential-based molecular dynamics (MD). DFT has been used to study the energetics and the electronic structure of the different dopant configurations, whereas MD is used to study thermal properties and oxygen diffusion at different temperatures.

2.1. Density functional theory calculations

DFT calculations were performed with the Vienna ab initio simulation package, VASP (version 5.4.1).^18–21 The projector-augmented wave method (PAW) was applied to describe the ion-electron interaction.²² Based on convergence tests, the plane wave energy cut-off for SmCoO₃ was set to 500 eV. Spin-polarized calculations were performed with the Perdew–Burke–Ernzerhof (PBE)^23,24 functional under electronic and ionic self-consistence, with convergence criteria of 10⁻⁵ eV and 10⁻³ eV Å⁻¹ respectively. The following valence electrons for the atomic species involved were taken into account: Sm (5s²5p⁶6s²), Sr (4s²4p⁶5s²), Ba (5s²5p⁶6s²), Ca (3s²3p⁶4s²), Co (4s²3d⁷), and O (2s²2p⁴). The tetrahedron method with Blöchl corrections for smearing²⁵ together with a 4 × 4 × 4 Γ-centered Monkhorst–Pack grid was applied.²⁶ Bader charges²⁷ were calculated using the Henkelman algorithm.²⁸ When non-stoichiometric bulk materials were modelled without charge compensation via oxygen vacancies, a homogeneous background charge was added to counteract the resulting charge. This makes the system overall charge-neutral and avoids a diverging Coulomb interaction.^29,30 To correct the electron self-interaction problem arising from the d-electrons of Co,^31–33 we have used the on-site Coulombic interaction (DFT+U), considering a Hubbard parameter of 3 eV using the Dudarev approximation.^34,35 Finally, the crystal structure used was the 2 × 2 × 2 Pm3̄m pseudo cubic cell (Fig. 1) as this has been shown to be sufficient to reproduce the properties of the system.³⁴ These materials have average cubic structures at high temperatures, with local Jahn–Teller distortions on the CoO-octahedra.³⁶

Fig. 1 Polyhedral representation of the 2 × 2 × 2 cubic SmCoO₃ perovskite structure. The pink atoms are Sm³⁺, whereas Co³⁺ are in the center of the blue octahedra, and O²⁻ are red. Sm³⁺ are also placed in the center of the edges and faces, but they are not shown here for the benefit of clarity. A²⁺ will be substituted at the Sm³⁺ lattice sites.

2.2. Molecular dynamics calculations

All MD simulations were performed using the DL_POLY 4.07 code,³⁷ with a 20 × 20 × 20 supercell simulated under NPT conditions, using a Nosé–Hoover thermostat.³⁸ All statistics were collected after a system equilibration of 13 ps, with an added 100 ps production phase. The time step for both phases was 0.5 fs. The Ewald summation was employed to account for the electrostatic interactions, whereas the Verlet algorithm evaluated the atomic motions. We have used the Buckingham potential, within the Born model for ionic solids,^39,40 to describe the short-range interactions between atoms.⁴¹ This potential is described in eqn (3), where r defines the distance between ions i and j, and A, ρ, and C are empirically fitted parameters.

In addition, we have considered the shell model to describe the electronic polarization of the atoms (α), where the ions are modeled as cores with a harmonic spring (k) connected to a massless shell with charge Y. The cores represent the nuclei and core electrons of the ions, whereas the shell represents the valence electrons. The relation is presented in eqn (4).

The interatomic potentials used in this work are taken from the Cherry et al. library, as they have previously been used to model ionic conductivity in perovskites.^42–45 The Sm³⁺–O²⁻ interatomic potential parameters had to be derived separately since they were not available for this potential model. To that end, we have used the GULP code^46–49 to fit them according to the structural data for cubic and orthorhombic SmCoO₃, as well as Sm₂O₃. The fitting results are presented in Table S1 of the ESI,† which show that the derived model perfectly replicates the structural properties of all different phases and materials considered. The full set of interatomic potentials used is presented in Table 1 below.

Table 1 Interatomic potential parameters used for simulation of Sm_1−xA_xCoO_3−x/2. Potential cutoff was set to 12 Å

Interaction	Short range interactions				Shell model
	A (eV)	ρ (Å)	C (eV Å^6)		Y (e)	k (eV Å^−2)	Ref.
Sm^3+–O^2−	1252.94	0.3590	0.00	Sm^3+	−0.250	173.90
Co^3+–O^2−	1329.82	0.3087	0.00	Co^3+	2.040	196.30	42
O^2−–O^2−	22764.30	0.1490	43.00	O^2−	−2.389	42.00	42
Ba^2+–O^2−	4818.42	0.3067	0.00	Ba^2+	1.831	34.05	44
Ca^2+–O^2−	1090.40	0.3437	0.00	Ca^2+	3.135	110.20	45
Sr^2+–O^2−	959.10	0.3721	0.00	Sr^2+	3.251	71.70	45

---

3. Results and discussion

3.1. Dopant configurations

To evaluate the position of the A²⁺ and the V_O, we have used the Site-Occupancy Disorder program (SOD).⁵⁰ SOD uses the symmetry operations of the non-doped bulk system to determine all the non-equivalent substitutions. At this stage, we should then use SOD again to determine the non-equivalent positions for the oxygen vacancy in the doped system. However, due to the high number of resulting configurations (around 250 configurations in total), we first determined the most stable dopant configuration in doped SmCoO₃ without V_O. We found three inequivalent configurations for x = 0.25, and six for x = 0.50 (ESI,† Fig. S1). Since SOD does not provide any information about the relative stability of the different bulks, we optimized all the different distributions using VASP, for all the different dopants: Ca²⁺, Sr²⁺, and Ba²⁺. After optimization, the most stable dopant distribution per dopant and concentration was used to determine the positions of the oxygen vacancies. We found that all systems with x = 0.25 had the same most stable dopant vacancy configuration, whereas Sr- and Ba-doped systems for x = 0.5 had different configurations from the Ca-doped system at the same concentration. The reason for this difference was attributed to size-effects, as Ca²⁺ is the smallest dopant. The final configurations for each cation and concentration are depicted in Fig. 2, and a full description of the energetics can be found in the ESI.†

Fig. 2 Dopant-vacancy configurations for (a) Sm_0.75A_0.25CoO_2.88, (b) Sm_0.5A_0.5CoO_2.75 (A = Ba, Sr), and (c) Sm_0.5Ca_0.5CoO_2.75. The pink spheres are Sm³⁺, whereas Co³⁺ are in the center of the blue octahedra, and O²⁻ are red. V_O is represented by white and red circles. In (a) grey spheres are Ca²⁺, Sr²⁺, or Ba²⁺, whereas in (b) green are Sr²⁺, or Ba²⁺, and in (c) light blue spheres are Ca²⁺.

3.2. Electronic and magnetic structures

Projected density of states (PDOS) for the oxidized systems (Sm_1−xA_xCoO₃) with x = 0.25 and 0.5 are collected in Fig. 3. All materials, with the exception of Sm_0.50Sr_0.50CoO₃, show the same common features. Upon doping, the systems become half-metallic, thus indicating that they all become electronic conductors. This conduction mainly occurs through the Co e_g orbitals, which are hybridized with the O 2p-orbitals. We modelled different magnetic structures for all the different systems, and in all cases the ferromagnetic ground state was the most stable. This magnetic structure is evidenced in the PDOS by the fact that α- and β-densities are not symmetrical. The change in the magnetic structure is attributed to the LS to IS transition observed in the cobalt centers and listed in Table 2, where Co³⁺ transitions from LS (μ_Co = 0 μ_B) to IS with μ_Co ≈ 2 μ_B. Cobalt charges (q_Co) do not vary markedly with increasing x, whereas q_O decreases by ≥0.04 e going from x = 0.25 to x = 0.5. This trend has also been observed for La_1−xSr_xFeO₃, and was suggested to indicate hole delocalization on the oxygen sub-lattice,¹⁶ although there is no clear evidence in the PDOS of the existence of this electronic hole near the Fermi level (E_F).

Fig. 3 PDOS for (a) Sm_0.75Ca_0.25CoO₃, (b) Sm_0.50Ca_0.50CoO₃, (c) Sm_0.75Sr_0.25CoO₃, (d) Sm_0.50Sr_0.50CoO₃, (e) Sm_0.75Ba_0.25CoO₃, and (f) Sm_0.50Ba_0.50CoO₃. Energies are referred to the Fermi level (E_F). All states below 0 eV represent occupied states, whereas above 0 eV, we find the virtual states. In addition, positive PDOS values are referred to α-spin occupations and negative values are referred to β-spin occupations.

Table 2 Bader charges (q_Co,O), and cobalt magnetic moments (μ_Co) for Sm_1−xA_xCoO_3−x/2. A is dopant, * denotes ion next to V_O (nearest neighbor). Sm, A, and O magnetic moments were found to be negligible in comparison to cobalt and are thus not included. Sm, and A Bader charges are included in the ESI. μ_Co,SmCoO3 = 0 μ_B, q_Co,SmCoO3 = 1.31 e, and q_O,SmCoO3 = −1.11 e

A		Sm0.75A0.25CoO3	Sm0.75A0.25CoO2.88	Sm0.5A0.5CoO3	Sm0.5A0.5CoO2.75
Ca	μ Co (μB)	1.8	0.13, 2.1	1.7	0.1, 2.0, 2.1
	μ Co* (μB)		1.8		1.9, 2.0, 2.6
	q Co (e)	1.36	1.21–1.32	1.36	1.23–1.32
	q Co* (e)		1.25		1.05–1.33
	q O (e)	−1.11 ± 0.02	−1.13 ± 0.1	−1.07 ± 0.02	−1.11 ± 0.07
Sr	μ Co (μB)	1.8	0.1, 2.1	0.01, 0.03, 0.04, 0.05	1.8, 2.0, 1.9, 2.2
	μ Co* (μB)		1.8		1.9, 2.1
	q Co (e)	1.38	1.23–1.32	1.33	1.32–1.37
	q Co* (e)		1.26		1.14–1.24
	q O (e)	−1.12 ± 0.02	−1.13 ± 0.06	−1.05 ± 0.02	−1.14 ± 0.06
Ba	μ Co (μB)	1.9	0.1, 2.1, 2.2	1.9	2.4, 3.1, 2.1, 0.4
	μ Co* (μB)		1.7		2.3, 2.0
	q Co (e)	1.40	1.22–1.30	1.40	1.28–1.48
	q Co* (e)		1.29		1.22–1.26
	q O (e)	−1.13 ± 0.03	−1.14 ± 0.1	−1.09 ± 0.04	−1.15 ± 0.1

---

For the specific case of Sm_0.50Sr_0.50CoO₃, it is still a conductor since the cobalt t_2g orbitals are found both above and below E_F, but its electronic structure reveals that the t_2g/e_g crystal field splitting still applies for this material despite the doping (Fig. 3d). Accordingly, the magnetic moment of cobalt atoms is about 0 and no significant change in the cobalt charge is observed.

For the V_O compensated systems, the PDOS are collected in Fig. 4. Now, without exception, all systems become half-metallic, thus indicating that they also become electronic conductors. The VBM is a mixture of O 2p, Co t_2g, and Co e_g, as opposed to the oxidized systems where Co e_g-states were mainly present at the CBM and above. Again, the conduction mainly occurs through the Co e_g-orbitals, which are hybridized with the O 2p-orbitals.

Fig. 4 PDOS for (a) Sm_0.75Ca_0.25CoO_2.88, (b) Sm_0.50Ca_0.50CoO_2.75, (c) Sm_0.75Sr_0.25CoO_2.88, (d) Sm_0.50Sr_0.50CoO_2.75, (e) Sm_0.75Ba_0.25CoO_2.88, and (f) Sm_0.50Ba_0.50CoO_2.75. Energies are referred to the Fermi level (E_F). All states below 0 eV represent occupied states, whereas above 0 eV, we find the virtual states. In addition, positive PDOS values are referred to α-spin occupations and negative values are referred to β-spin occupations.

The presence of V_O induces different magnetic moments in the cobalt centers. Those who are next to the V_O keep the same magnetic moment observed in the oxidized systems, around 1.8 μ_B, whereas the rest of the cobalt atoms show different values, ranging from 0 to 2.2 μ_B. All magnetic moments show spin in the same direction, indicating that the doped systems with oxygen vacancies become ferromagnetic (Table 2). Sm_0.50Ba_0.50CoO_2.75 has the highest μ_Co, indicating that Co³⁺ is in the IS/HS spin state. This mixing of μ_Co is in agreement with experimental results on doped SmCoO₃.⁸

The presence of V_O in the lattice gives different cobalt spin states within the system, but not cobalt mixed valence, which has also been noted experimentally for x > 0.2.⁵¹ This finding would imply that the charge compensation is fully completed through V_O compensation, without remaining negative charge in V_O. A similar result was found for V_O formation in the undoped SmCoO₃.³⁴ However, the subtle difference in q_Co between the oxidized and V_O compensated systems is worth noting. The difference in total q_Co in the two different charge compensation schemes for the same A and x is, however, larger. Δ_Ca,0.25 = 0.70 e, Δ_Ca,0.50 = 0.91 e, Δ_Sr,0.25 = 0.78 e, Δ_Sr,0.50 = 0.33 e, Δ_Ba,0.25 = 0.92 e, and Δ_Ba,0.50 = 0.65 e. This could possibly indicate the presence of an electron hole on the cobalt atoms, but we would then expect charges closer to 2 e. Therefore, we consider that due to how the Bader charges are partitioned, and that Δ_A,x do not sum to a whole electron, no certain conclusion can be drawn.

Since the analysis of the density of states, and the variation of the cobalt magnetic moments and Bader charges do not provide clear evidence of the actual existence of an electron hole to compensate the charge of the system, we have examined their possible presence by calculating the charge density difference (Δρ) according to eqn (5).

Δρ = ρ_{Sm1−xAxCoO3−x/2} + nρ₀ − ρ_{Sm1−xAxCoO3} (5)

where ρ is the charge density and n the number of V_O in lattice. If there is an electronic hole, it should be represented as a depletion in the charge density plot. According to Fig. 5, however, Δρ is mainly centered on the oxygen ions, indicating that the variation in charge density is located on the oxygen sub-lattice. Although we can observe charge depletion close to the cobalt centers, this is not enough to fully verify the presence of electron holes on the cobalt atoms. Furthermore, rather than an electron hole, the charge density difference illustrates a rearrangement of charges rather than a proper hole. We also considered hole compensation in terms of the formation of O⁻, by comparing the lattice relaxation around the oxygen nearest to the dopants. This formation energy can be obtained by calculating the difference in total energy between the fully relaxed localized polaron state, and the hole state in the perfect undistorted lattice, i.e. unrelaxed system.^52,53 From these calculations, we found that the hole polaron formation energy is 0.62 eV for x_Ca = 0.25; 0.95 eV for x_Ca = 0.50; 0.92 eV for x_Sr = 0.25; 0.78 eV for x_Sr = 0.50; 1.06 eV for x_Ba = 0.25; and 0.60 eV for x_Ba = 0.50. These formation energies are quite high, although could become favorable under certain conditions.

Fig. 5 Charge density difference for (a) x_Ca = 0.25, (b) x_Ca = 0.5, (c) x_Sr = 0.25, (d) x_Sr = 0.5, (e) x_Ba = 0.25, and (f) x_Ba = 0.5. Isosurface value is 0.5. Yellow isosurface is positive Δρ (gain of charge), and blue is negative Δρ (depletion of charge).

The question still remains though, why A²⁺-doped SmCoO₃ becomes an electronic conductor. The key to understand this behavior can be seen in Sm_0.50Sr_0.50CoO₃, which is almost diamagnetic and still shows e_g/t_2g splitting. Examining the geometric structure of Sm_0.50Sr_0.50CoO₃, we realize that there are no distortions in the Co–O bond lengths (see full list in Table S4 in ESI†), and thus no CoO₆-octahedral distortions. Therefore, the crystal field splitting is still intact, which explains why Sm_0.50Sr_0.50CoO₃ has not become at least half metallic. From this finding, we can conclude that any doping that induces structural distortion in SmCoO₃, will alter the electronic structure, i.e. will break the crystal field splitting, and the Co d-orbitals will be closer in energy, as well as more hybridized with O 2p-orbitals, turning the system into an electrical conductor. To prove this, hypothesis we calculated the 0 K electronic conductivity (σ_e) for all Sm_1−xA_xCoO₃, and plotted them as a function of Co–O bond length variance. The variance was set as the difference between the shortest and the longest Co–O distance, according to Table S4 in ESI.† The larger the spread in Co–O bond lengths, the higher the electronic conductivity. The trend is presented in Fig. 6. Although it is not a perfect linear dependence, it is indicative enough to conclude that conductivity in Sm-site doped SmCoO₃ is a direct consequence of the distortion induced by doping the system, rather than being a consequence of any electron hole.

Fig. 6 Electronic conduction plotted against Co–O bond length variance within the systems.

3.3. Oxygen vacancy formation energy

To calculate (Table 3), we have employed the method developed by Zhang and Northrup (eqn (6)),^54,55where E_Defective is the total energy of the defective system, E_Perfect is the total energy of the non-defective system; n_O is the number of removed or added oxygen ions from the bulk, and μ_O is the oxygen chemical potential.^{16,34,56–60} For V_O we consider an O-rich situation, which for SmCoO₃ equals to μ_O = ½E_O2 − 0.5, as it was demonstrated in the phase diagram for this perovskite.³⁴ This approach has been successfully employed to perovskite systems in other computational studies.¹⁵ It is worth noting that we assume that the introduction of dopants in the system will not affect the chemical potential for oxygen. Furthermore, DFT overestimates the oxygen binding energy, and by how much is dependent on several computational parameters.⁶¹ However, inclusion of the oxygen overpotential correction in eqn (6), based on corrections used in previous studies,¹⁷ does not influence the results presented here, i.e. all oxygen vacancy formation energies are negative, and much lower in energy than the hole compensation scheme^34,61 Furthermore, we have neglected thermal, vibrational, and entropic contributions, as these are commonly known to be smaller than the typical DFT error.^36,62 It should also be noted that in real materials these properties are dependent on both temperature and oxygen partial pressure, which could directly influence these values. The dependence on temperature and partial oxygen pressure of the oxygen vacancy concentration and oxygen vacancy formation energy, which would affect the abundance of holes in these materials (eqn (1) and (2)), have been widely discussed elsewhere,^{36,54,58,61,62} and is thus only mentioned here for the sake of context.

Table 3 Oxygen vacancy formation energy for the most stable configuration. Other V_O-configurations are included in ESI

System	, x = 0.25 (eV)	, x = 0.50 (eV)
CaxSm1−xCoO3−x/2	−1.39	−0.44
SrxSm1−xCoO3−x/2	−0.95	−0.49
BaxSm1−xCoO3−x/2	−0.53	−1.11

---

V_O formation in these doped materials is spontaneous, i.e. negative , as expected considering the oxygen vacancy charge compensation scheme. Moreover, as mentioned for the hole compensation, dopant-size effects appear to exist when the systems obtain more ACoO₃ character at the higher dopant concentrations, leading to opposing trends in at the higher concentration. This trend has previously been observed for AMnO₃, with decreasing going from CaMnO₃, via SrMnO₃, and to BaMnO₃.⁶³ However, as all are negative, the effect would not be noticeably different between the different materials. Furthermore, the negative does indicate that oxygen could spontaneously leave the lattice, which would eventually lead to breakdown of the material. Material breakdown is known in the literature to occur at high temperatures, which are, however, outside the IT-SOFC range.^64,65 It is also worth noting that the oxygen vacancy formation energies (corresponding to eqn (1)) are much lower in all cases than the hole formation energies (corresponding to eqn (2)), which would indicate that oxygen vacancy formation is the most probable compensation scheme under equilibrium conditions.

Compared to SmCoO₃, which had an , we can conclude that doping clearly favors the formation of spontaneous vacancies. It is also worth mentioning that at high temperatures, SmCoO₃ has a C-type antiferromagnetic structure.⁶⁸ Thus, we also calculated for the high-temperature structure, and found this to be 1.43 eV. Furthermore, for comparison, we calculated for the traditional high-temperature SOFC cathode material La_0.75Sr_0.25MnO_3−d (LSM), which is 1.82 eV. We thus see that V_O formation in Sm_1−xA_xCoO_3−x/2 is energetically more favorable than in LSM. Thus, all of these systems show excellent potential as ionic conductors, but these results are not conclusive as to which x or A²⁺ leads to the highest oxygen conduction. Next, we have therefore used MD to evaluate the oxygen diffusion and ionic conduction. Please note that all the results below are for the systems in the V_O compensated scheme.

3.4. Association energy

Another important property to assess is the dopant-V_O-cluster interaction energy (E^int), which is calculated (eqn (7)) as the difference in E_f of the defect cluster (E^pair_f), and the sum of E_f for the two individual defects (E^{isolated defects}_f).⁵⁸

E^int = E^pair_f − ∑E^{isolated defects}_f (7)

Negative values (e.g. −0.67 eV for Ca-doping) indicate an energetic benefit for the defects to be associated, whereas positive values (e.g. 0.06 eV for Sr-doping and 0.20 eV for Ba-doping) indicate the tendency for the defects to dissociate. Thus, trapping of V_O may occur near a Ca-dopant site at low temperatures, although at operating temperatures, either intermediate or high, this barrier will be easy to overcome and should not noticeably influence oxygen migration.

3.5. Oxygen migration and conduction

Oxygen diffusion in perovskite oxides undergoes a vacancy hopping mechanism and is temperature dependent. Diffusion can be calculated using Fick's first law, assuming steady state diffusion. For steady state diffusion, the diffusion flux (J), which is the number of atoms diffusing per unit area and time, along a given direction (x, y, or z) is proportional to the concentration (C) gradient (eqn (8)).where D_i is the self-diffusion coefficient for species i, defined as (eqn (9))D₀ is the temperature-independent pre-exponential, E_a is the diffusion activation energy, k_B Boltzmann's constant, and T temperature. E_a is thus obtained from the gradient of ln(D_i) against 1/T Arrhenius plots. To obtain D_i, the mean square displacement (MSD) is calculated (using the Einstein relation) for the different species, and it is then obtained from the gradient of the MSD plot through eqn (10).

〈r²(t)〉 = 6D_it + B_i (10)

where 〈r²(t)〉 is the time-dependent MSD, t is time, and B_i is the thermal factor associated with atomic vibrations.

According to the MSD representation depicted in Fig. S1 of the ESI,† cation diffusion is limited and reaches equilibrium shortly after equilibration, whereas the oxygen MSD clearly increases with time (Fig. S2 in ESI†). D_O for oxygen was calculated according to eqn (10), and showed that in general the highest oxygen conductivity is found at temperatures higher than 1000 K, and for Ca_xSm_1−xCoO_3−x/2 at 1500 K. Moreover, all materials show higher oxygen diffusion (even at low temperatures) than the traditional cathode material La_0.5Sr_0.5MnO_3−d, which has a experimentally measured D_O of 3 × 10⁻¹² cm² s⁻¹ at 900 °C (1173 K),⁶⁹ as compared to the calculated value of 2.87 × 10⁻⁷ cm² s⁻¹ for Sm_0.5Ca_0.5CoO_2.75 at 1200 K (Table S5, ESI†).

The highest oxygen diffusion is reached when the system is doped with Ca²⁺, which can be explained by the fact that Ca²⁺ is also the system with the smallest of the tested dopants. Relating diffusion to activation energy (Fig. 7 and Table 4), E_a decreases with increasing x, indicating that higher dopant concentrations favor oxygen migration, as has previously been reported for Sr-doping of LaMnO₃ and LaCoO₃.⁴¹

Fig. 7 (a) Arrhenius plots for oxygen diffusion, and (b) ionic conductivity for Ca_xSm_1−xCoO_3−d, Ba_xSm_1−xCoO_3−d and Sr_xSm_1−xCoO_3−d. Temperature range is 800–1500 K.

Table 4 Oxygen diffusion activation energies (eV) for temperature range 800–1500 K

System	E a, x = 0.25 (eV)	E a, x = 0.50 (eV)
BaxSm1−xCoO3−d	0.42	0.32
CaxSm1−xCoO3−d	0.45	0.36
SrxSm1−xCoO3−d	0.42	0.35
SmCoO3	0.93

---

Ionic conduction (σ_O) is challenging to measure experimentally since both electronic and ionic conductivities are normally hard to differentiate.¹⁰ Thus, theoretical calculations have proved to be very useful in directly comparing σ_O between different systems. σ_O is calculated from D_O through the Nernst–Einstein equation¹⁰ (eqn (11)).

where C_O is the oxygen concentration, and q the charge of oxygen. The ionic conductivities of the different systems are plotted in Fig. 7b, which are comparable to experimental studies conducted by Yeh et al. on Sm_0.5Sr_0.5CoO_2.75, which gave log(σ_OT) between 0 and 2 for a temperature range of 700–1000 K (please note these values have been read from a graph).¹⁰ Furthermore, the ionic conductivities of these materials are close to those of SOFC electrolytes GDC and YSZ, which are excellent oxygen conductors.¹⁰ Thus, all systems investigated here have good oxygen conductivity for SOFC purposes. Specifically, the systems showing highest ionic conductivity at IT-SOFC temperatures are Sm_0.5Sr_0.5CoO_2.75, Sm_0.5Ca_0.5CoO_2.75 and Sm_0.75Ca_0.25CoO_2.83.

3.6. Thermal expansion coefficient

In addition to ionic conductivity, another vital property of SOFC cathodes is the thermal expansion coefficient (TEC). The cathode requires a TEC that is compatible with the electrolyte, as incompatibility will lead to SOFC cell breakdown. Common SOFC electrolytes include YSZ, LSGM, and GDC, which have TEC of 10.0–12.5 × 10⁻⁶ K⁻¹.⁷⁰ TEC for SmCoO₃ (17.7 × 10⁻⁶ K⁻¹) itself is incompatible with these electrolytes, and thus to reduce TEC, doping is vital.⁶⁷ DFT+U results showed that the lattice volume increases with increasing x, thus inducing negative chemical pressure, as all dopants have a higher ionic radius than Sm³⁺; Sm³⁺(0.96 Å) < Ca²⁺(1.00 Å) < Sr²⁺(1.18 Å) < Ba²⁺(1.35 Å).⁷¹ The largest volume increase for all x is seen for Ba²⁺, and the least for Ca²⁺, which behavior has previously been reported in molecular dynamics studies of Sr-doped SmCoO₃.⁷² The volume expansion in relation to temperature was then calculated with MD, and the linear thermal expansion coefficient extracted (Table 5).

Table 5 Linear thermal expansion coefficient. Orthorhombic SmCoO₃ has an experimental linear thermal expansion coefficient of 21.7 × 10⁻⁶ K^−1 66

System	TEC, x = 0.25 (10^−6 K^−1)	TEC, x = 0.50 (10^−6 K^−1)
CaxSm1−xCoO3−d	23.3	24.2
SrxSm1−xCoO3−d	22.7	30.9
BaxSm1−xCoO3−d	23.7	30.0

---

System	TEC, x = 0 (10^−6 K^−1)
SmCoO3	17.1

---

The TECs are in alignment with previous results for cobalt-based SOFC cathodes, which normally have higher TEC than SOFC electrolytes.⁷⁰ To overcome this incompatibility between these cathodes and traditional electrolytes, Co-site doping with transition metals might be necessary. Such doping has been shown to decrease cathode TEC in for example La_1−xSr_xCoO_3−x/2.⁷³ However, this is outside the scope of this paper.

Comparing all the above calculated properties of Sm_1−xA_xCoO_3−x/2, it can be seen that Ca-doping gives the highest D_O, among the lowest E_a, highest ionic conduction, and more compatible TEC at intermediate temperatures (800–1000 K). Sr-Doping, however, gives higher ionic conduction at lower temperature (x = 0.50), but due to this composition a much higher TEC. Thus, Ca-doping might be more practical for IT-SOFC cathode applications.

4. Conclusion

The electronic properties and ionic conduction of potential IT-SOFC cathode materials have been modeled using DFT+U and MD calculations. Evaluation of two charge compensation schemes suggests that oxygen vacancy compensation is the more likely for these materials. Furthermore, the change in electronic structure upon doping has been linked directly to distortions in the Co–O bonds. Structural distortions break the octahedral crystal field splitting, which provokes a spin transition in Co that leads to a change in the electronic structure, making the doped SmCoO₃ a conductor.

However, we have found that the ionic conductivity in doped SmCoO₃ is two times higher at intermediate temperatures than in LSM, with Ca²⁺, x = 0.25, reporting the best results, in agreement with experimental findings. Finally, thermal expansion coefficients were also calculated, and were found to be in the same range as lanthanum-based cathode materials, but still too high for perfect compatibility with the traditional electrolyte materials. One way to overcome this deficiency could be the introduction of additional Co-site dopants, which is the topic of future studies.

Acknowledgements

The authors acknowledge the Engineering and Physical Sciences Research Council (EPSRC) for financial support (Grants references EP/K016288/1 and EP/K001329/1). We also acknowledge the use of the UCL Legion High Performance Computing Facility (Legion@UCL), as well as the use of the UCL Grace High Performance Computing Facility (Grace@UCL), and associated support services, in the completion of this work. Finally, via our membership of the UK's HPC Materials Chemistry Consortium, which is funded by EPSRC (EP/L000202), this work made use of the facilities of ARCHER, the UK's national high-performance computing service, which is funded by the Office of Science and Technology through EPSRC's High End Computing Programme. DL_POLY_4 is a molecular dynamics simulation package written by I. T. Todorov and W. Smith, and has been obtained from STFC's Daresbury Laboratory via the website http://www.ccp5.ac.uk/DL_POLY. N. H. d. L. thanks the Royal Society for an Industry Fellowship. E. O. gratefully acknowledges EPSRC funding from the Centre for Doctoral Training in Molecular Modelling and Materials Science (EP/G036675/1).

References

1. S. M. Haile, Acta Mater., 2003, 51, 5981–6000.
2. M. M. Kuklja, E. A. Kotomin, R. Merkle, Y. A. Mastrikov and J. Maier, Phys. Chem. Chem. Phys., 2013, 15, 5443–5471.
3. B. C. Steele and A. Heinzel, Nature, 2001, 414, 345–352.
4. A. Aguadero, L. Fawcett, S. Taub, R. Woolley, K.-T. Wu, N. Xu, J. A. Kilner and S. J. Skinner, J. Mater. Sci., 2012, 47, 3925–3948.
5. Y. Choi, D. S. Mebane, M. C. Lin and M. Liu, Chem. Mater., 2007, 19, 1690–1699.
6. H. Fukunaga, M. Koyama, N. Takahashi, C. Wen and K. Yamada, Solid State Ionics, 2000, 132, 279–285.
7. X. Jiang, Q. Xu, Y. Shi, X. Li, W. Zhou, H. Xu and Q. Zhang, Int. J. Hydrogen Energy, 2014, 39, 10817–10823.
8. T. N. Vasil’chikova, T. G. Kuz’mova, A. A. Kamenev, A. R. Kaul’ and A. N. Vasil’ev, JETP Lett., 2013, 97, 34–37.
9. I. Fullarton, J. Jacobs and H. Van Benthem, Ionics, 1995, 1, 51–58.
10. T. C. Yeh, J. L. Routbort and T. O. Mason, Solid State Ionics, 2013, 232, 138–143.
11. C. Xia, W. Rauch, F. Chen and M. Liu, Solid State Ionics, 2002, 149, 11–19.
12. B. Wei, K. Chen, C. C. Wang, Z. Lü and S. P. Jiang, Electrochim. Acta, 2015, 174, 327–331.
13. S. Pawar and K. Shinde, J. Mater., 2013, 2013, 6.
14. F. A. Kröger and H. J. Vink, Solid State Phys., 1956, 3, 307–435.
15. A. M. Deml, V. Stevanović, C. L. Muhich, C. B. Musgrave and R. O’Hayre, Energy Environ. Sci., 2014, 7, 1996.
16. A. M. Ritzmann and A. B. Muñoz-García, Chem. Mater., 2013, 25, 3011–3019.
17. A. M. Ritzmann, J. M. Dieterich and E. A. Carter, Phys. Chem. Chem. Phys., 2016, 18, 12260–12269.
18. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 47, 558.
19. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 49, 14251.
20. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15–50.
21. G. Kresse and J. Furthmüller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169–11186.
22. P. E. Blöchl, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 17953.
23. J. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
24. J. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1997, 78, 1396.
25. P. E. Blöchl, O. Jepsen and O. K. Andersen, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 49, 16223–16233.
26. H. J. Monkhorst and J. D. Pack, Phys. Rev. B: Solid State, 1976, 13, 5188–5192.
27. R. F. W. Bader, Atoms in Molecules: A Quantum Theory, Oxford University Press, Oxford, 1990.
28. G. Henkelman, A. Arnaldsson and H. Jónsson, Comput. Mater. Sci., 2006, 36, 354–360.
29. G. Pacchioni, J. Chem. Phys., 2008, 128, 182505.
30. H.-P. Komsa, T. T. Rantala and A. Pasquarello, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 86, 45112.
31. T. Geng, Z. Han and S. Zhuang, Phys. B, 2010, 405, 3714–3716.
32. P. Ravindran, A. Kjekshus, H. Fjellvåg, A. Delin and O. Eriksson, Phys. Rev. B: Condens. Matter Mater. Phys., 2002, 65, 64445.
33. P. Ravindran, R. Vidya, H. Fjellvåg and A. Kjekshus, J. Cryst. Growth, 2004, 268, 554–559.
34. E. Olsson, X. Aparicio-Anglès and N. H. De Leeuw, J. Chem. Phys., 2016, 145, 14703.
35. S. Dudarev, G. Botton and S. Savrasov, Phys. Rev. B: Condens. Matter Mater. Phys., 1998, 57, 1505–1509.
36. S. Zhang, N. Han and X. Tan, RSC Adv., 2015, 5, 760.
37. I. T. Todorov, W. Smith, K. Trachenko and M. T. Dove, J. Mater. Chem., 2006, 16, 1911–1918.
38. W. G. Hoover, Phys. Rev. A: At., Mol., Opt. Phys., 1985, 31, 1695–1697.
39. M. Born and J. Mayer, Phys. B, 1932, 75, 1–18.
40. J. E. Mayer, J. Chem. Phys., 1933, 1, 270–279.
41. M. S. Islam, M. Cherry and C. R. A. Catlow, J. Solid State Chem., 1996, 124, 230–237.
42. M. Cherry, M. S. Islam and C. R. A. Catlow, J. Solid State Chem., 1995, 118, 125–132.
43. M. S. D. Read, M. S. Islam, G. W. Watson, F. King and F. E. Hancock, J. Mater. Chem., 2000, 10, 2298–2305.
44. T. S. Bush, J. D. Gale, C. R. A. Catlow and P. D. Battle, J. Mater. Chem., 1994, 4, 831.
45. G. V. Lewis and C. R. A. Catlow, J. Phys. C: Solid State Phys., 1985, 18, 1149–1161.
46. J. D. Gale, J. Chem. Soc., Faraday Trans., 1997, 93, 629–637.
47. J. D. Gale and A. L. Rohl, Mol. Simul., 2003, 29, 291–341.
48. J. D. Gale, Z. Kristallogr., 2005, 220, 552–554.
49. J. D. Gale, Philos. Mag. B, 1996, 73, 3–19.
50. R. Grau-Crespo, S. Hamad, C. R. A. Catlow and N. H. De Leeuw, J. Phys.: Condens. Matter, 2007, 19, 256201.
51. Y. Orikasa, T. Ina, T. Nakao, A. Mineshige, K. Amezawa, M. Oishi, H. Arai, Z. Ogumi and Y. Uchimoto, J. Phys. Chem. C, 2011, 115, 16433–16438.
52. A. L. Shluger, K. P. McKenna, P. V. Sushko, D. M. Ramo and A. V. Kimmel, Modell. Simul. Mater. Sci. Eng., 2009, 17, 84004.
53. E. A. Kotomin, R. I. Eglitis, A. V. Postnikov, G. Borstel and N. E. Christensen, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 60, 11–15.
54. S. T. Murphy and N. D. M. Hine, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 87, 1–18.
55. S. Zhang and J. Northrup, Phys. Rev. Lett., 1991, 67, 2339–2342.
56. H. Raebiger, S. Lany and A. Zunger, Phys. Rev. B: Condens. Matter Mater. Phys., 2007, 76, 45209.
57. T. Tanaka, K. Matsunaga, Y. Ikuhara and T. Yamamoto, Phys. Rev. B: Condens. Matter Mater. Phys., 2003, 68, 205213.
58. P. G. Sundell, M. E. Björketun and G. Wahnström, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 73, 104112.
59. Z. Alahmed and H. Fu, Phys. Rev. B: Condens. Matter Mater. Phys., 2007, 76, 224101.
60. E. Heifets, S. Piskunov, E. A. Kotomin, Y. F. Zhukovskii and D. E. Ellis, Phys. Rev. B: Condens. Matter Mater. Phys., 2007, 75, 115417.
61. Y.-L. Lee, J. Kleis, J. Rossmeisl and D. Morgan, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 80, 224101.
62. J. Buckeridge, D. O. Scanlon, A. Walsh and C. R. A. Catlow, Comput. Phys. Commun., 2014, 185, 330–338.
63. A. M. Ritzmann, M. Pavone, A. B. Muñoz-García, J. A. Keith and E. A. Carter, J. Mater. Chem. A, 2014, 2, 8060–8074.
64. C. Persson, Y. J. Zhao, S. Lany and A. Zunger, Phys. Rev. B: Condens. Matter Mater. Phys., 2005, 72, 1–14.
65. A. Marthinsen, C. Faber, U. Aschauer, N. A. Spaldin and S. M. Selbach, MRS Commun., 2016, 1–26.
66. R. Liu, D. Xu, S. Li, Z. Lü, Y. Xue, D. Wang and W. Su, Front. Chem. China, 2006, 1, 398–401.
67. S. W. Baek, J. H. Kim and J. Bae, Solid State Ionics, 2008, 179, 1570–1574.
68. E. Olsson, X. Aparicio-Anglès and N. H. de Leeuw, J. Chem. Phys., 2016, 145, 224704.
69. S. P. Jiang, J. Mater. Sci., 2008, 43, 6799–6833.
70. J. A. Kilner and M. Burriel, Annu. Rev. Mater. Res., 2014, 44, 365–393.
71. R. D. Shannon, Acta Crystallogr., Sect. A: Cryst. Phys., Diffr., Theor. Gen. Crystallogr., 1976, 32, 751–767.
72. M. A. Farhan and M. J. Akhtar, J. Phys.: Condens. Matter, 2010, 22, 75402.
73. J. H. Kim, S.-W. Baek, C. Lee, K. Park and J. Bae, Solid State Ionics, 2008, 179, 1490–1496.

---

Footnote

† Electronic supplementary information (ESI) available: Calculated and experimental lattice parameters for cubic and orthorhombic SmCoO₃, as well as Sm₂O₃ using interatomic potentials (Table S1), a description of the Goldschmidt factor and its calculated values (Table S2), the relative energies of the different dopant configurations (Table S3), lattice parameters, and interatomic distances in the dopant systems for the most stable configurations (Table S4), and a discussion of the solution energy (Table S5). Furthermore, graphical representations of all oxygen vacancy-dopant configurations and their oxygen vacancy formation energies are included in Fig. S1–S3. In Tables S6 and S7, Bader charges have been collected, and mean square displacement plots are presented in Fig. S4 and S5. Table S8 contains all diffusion coefficients used to calculate ionic conductivity and activation energies. See DOI: 10.1039/c7cp01555k

---