Article pubs.acs.org/JPCC
Improving Oxygen Transport in Perovskite-Type LaGaO3 Solid Electrolyte through Strain Cristina Tealdi* and Piercarlo Mustarelli Department of Chemistry and INSTM, University of Pavia, Viale Taramelli 16, 27100 Pavia, Italy S Supporting Information *
ABSTRACT: Lattice strain is a promising possibility to improve materials performance in view of their application in thin-film devices. In particular, defect and transport properties in ionic conductors may be tailored through strain effects, since defect formation energy and migration barriers are correlated to structural parameters which, in turn, are influenced by strain-induced deformations. In this computational study we predicted that oxide-ion diffusion in perovskite-type lanthanum gallate can be improved through application of tensile strain. The structural deformations required to accommodate tensile lattice strain in the perovskite system are shown to result in a preferential localization of the oxygen vacancies in the equatorial plane of the GaO6 octahedra, while oxide-ion diffusion becomes anisotropic.
1. INTRODUCTION Ionic conduction is at the basis of relevant “clean energy” technologies, one of the most important being indeed that of solid oxide fuel cells (SOFCs). In such devices, stringent requirements apply on the ionic conductivity of the electrolyte layer, which at the operating temperature of interest should be at least 10−2 S cm−1 in order for the material to be technologically competitive.1 While the search for new compositions and structures capable of delivering high ionic conductivity at intermediate temperature (500 < T < 700 °C) with reduced costs and environment impact is extremely active,2−8 the most used electrolyte materials to date are still yttrium-stabilized zirconia (YSZ), gadolinium-doped ceria (CGO), and perovskite-type lanthanum gallate (LSGM).1,9,10 Whereas the role of both the nature and the amount of dopants was thoroughly investigated in the past, nanostructuring and reduction of the electrolyte thickness are currently considered as the most interesting routes to increase material performance.11 In particular, the report of a colossal conductivity increase in epitaxial heterostructures formed by thin layers of YSZ and strontium titanate (STO)12 prompted many investigations toward substrate-induced strain effects on the transport properties of functional oxides from both the computational and the experimental points of view.13−20 In this regard, from the computational side, the main focus was put on the effects of tensile and compressive strain on the transport and defect properties of oxide-ion conductive materials. In particular, the defect and transport properties of YSZ13−15 and CGO16−18 under mechanical stress were studied by means of static lattice, molecular dynamics, and/or ab initio calculations, focusing attention on aspects such as oxygen diffusion, dopant-vacancy association energies, and migration barriers. These studies generally showed that the application of © XXXX American Chemical Society
strain, especially tensile strain, has a positive effect on the ionic conductivity of fluorite-type oxide-ion conductors. They also pointed out that the extraordinary increase of 8 orders of magnitude reported for the epitaxial YSZ/STO heterostructures is difficult to be explained based on computational results as well as on structural and thermodynamic considerations.17 On the other hand, since such studies did not explicitly include interfacial phenomena but focused on the strain effects in the bulk, comparison with experimental data on complex heterostructures requires some attention, bearing in mind the model limitations compared to the complex behavior of real systems. However, these studies showed the importance of analyzing strain effects, e.g., strain and orientation due to epitaxial growth, as a mean to tailor defects in the chemistry and ionic mobility in thin-film devices. The present work offers a wide computational investigation of the structural, defect, and transport properties of perovskitetype LSGM, thus complementing the survey of strain effects on such properties related to well-established electrolyte materials for SOFCs already present in the literature for fluorite-type materials (e.g., YSZ and CGO). In particular, the focus is on the effects of biaxial tensile or compressive strain in the [100] and [010] directions of orthorhombic LaGaO3. Static lattice and molecular dynamics techniques are used complementary. Despite LSGM being a widely studied electrolyte material, this work represents, to the best of our knowledge, the first application of classical molecular dynamics to the investigation of oxygen vacancies diffusion in this system. In fact, a previous Received: September 17, 2014 Revised: October 24, 2014
A
dx.doi.org/10.1021/jp509413w | J. Phys. Chem. C XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry C
Article
Figure 1. Evolution of the pseudocubic unit cell parameters and volume as a function of biaxial applied strain for orthorhombic LaGaO3.
minimized structure at constant pressure in condition of zero applied strain was used to generate biaxial compressive and tensile strains between −5% and +5% in the ab plane. Energy minimization of the structure with respect to the c cell parameter and all atomic coordinates was then allowed. The strained energy-minimized structures were then used as the starting point for subsequent defect calculations. For MD simulations, a doped composition was considered: La0.8Sr0.2GaO2.9. In this case, dopants and oxygen vacancies were randomly introduced in a 6 × 6 × 6 supercell containing more than 4000 atoms by respecting charge neutrality. As for the pure system, the energy-minimized structure at constant pressure in condition of zero applied strain was used to generate biaxial compressive and tensile strains between −4% and +5% in the ab plane. As for the undoped system, energy minimization of the structure with respect to the c cell parameter and all atomic coordinates was then allowed for each applied strain value based on SL simulations. The strained energy-minimized structures were then used as the starting point for molecular dynamic simulations. The short-range parameters used for SL simulations were transferred to the MD simulations that were performed with the Dl_Poly code.26 The shell model was not implemented in the MD simulations. The simulation temperature was 1873 K. Such a high simulation temperature was chosen to allow rapid and efficient thermalization of the highly defective solid; it also allowed sufficient oxygen diffusion to occur during the time span covered by the MD simulation, thereby improving the statistics of the calculations. The unstrained system was first equilibrated at the desired temperature under a constant pressure of 1 atm for 100 000 time steps (NPT ensemble) with a time step of 1 fs by scaling the temperature every 10 steps. After thermalization, the simulation box was equilibrated for at least further 30 000 steps with box lengths held fixed at their average values obtained during the NPT stage. The main simulation run of 300 000 time steps was then performed in the NVT ensemble (Nose−Hoover thermostat) to give a total simulation time of 300 ps, with data recorded every 100 time steps for postsimulation analysis. For the strained systems, the energy-minimized cell parameters, obtained from SL simulations for each applied strain, were expanded according to the thermal expansion coefficient derived for the same composition in the absence of applied strain. Simulations in the NVT ensemble were then performed as described above.
molecular dynamics work dealt with diffusion of cationic species in LSGM.21 As for previous computational studies dealing with the effect of strain on the transport properties of SOFCs electrolyte materials,12−17 interfaces phenomena are not included explicitly, and only the bulk effects of the lattice strain are taken into account. Therefore, we stress again that also our results may offer a partial explanation of the complex behavior of substrateinduced strain effects in thin films. However, they are important in elucidating the atomic-scale phenomena influencing structural, defect, and transport properties of functional oxides in thin-film devices. We are confident that the predictive nature of the present work will shed light on oxygen transport in strained complex structures and stimulate development of novel functional layers based on LaGaO3.
2. METHODS This study was performed through the use of two simulation techniques, namely, static lattice (SL) and molecular dynamics (MD). In brief, SL simulations involve calculation of the potential energy of the system in terms of atomic coordinates, with ionic interactions described as ionic pairwise potentials. The interactions between ions are partitioned into long-range Coulombic forces and short-range forces accounting for electron cloud overlap (Pauli repulsion) and van der Waals interactions. The short-range pair potential interactions were modeled with a Buckingham potential, whose parameters were obtained by previous computational work on this system.22 Electronic polarizability of the ions was taken into account using the shell model developed by Dick and Overhauser.23 The energy of the charged defect relative to the perfect, energyminimized lattice was calculated using the Mott−Littleton approach according to which the lattice is partitioned into inner- and outer-spherical regions centered on the defect or defect cluster. Ions in the inner region were relaxed explicitly, whereas the remainder of the crystal was treated by more approximate quasi-continuum methods. SL simulations were performed with the use of the General Utility Lattice Program (GULP).24 Structural data for the LaGaO3 system in orthorhombic symmetry (space group Pbnm) were used as the starting point for the calculations.25 In order to introduce biaxial strain in the ab plane and to evaluate its effect on the structural and transport properties of the LaGaO3 system, the following procedure was adopted for the pure system. The energyB
dx.doi.org/10.1021/jp509413w | J. Phys. Chem. C XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry C
Article
for this system22 support the feasibility of applying the present set of potential parameters to the computational investigation of perovskite-type LaGaO3. Careful analysis of the structural evolution of selected parameters as a function of the applied strain can give information on how the perovskite structure is able to accommodate compressive and tensile strain through tilting and deformation of the BO6 octahedra. Such structural information is relevant to the comprehension and optimization of transport properties. Figure 3 shows the energy-minimized structure of LaGaO3 for −5%, 0%, and +5% applied strain in the ab plane. Visual
The diffusion coefficient, D, of a species was calculated from the slope of the mean square displacement (MSD) vs time plot using the Einstein relation
d(MSD) = B + 6Dt dt
(1)
where B is twice the square of the mean vibrational amplitude of the species. The nMoldyn software was used to analyze the MD data.27
3. RESULTS AND DISCUSSION 3.1. Structural Properties. The first information that can be extracted from analysis of the SL results is concerned with structural parameter changes in the pure undoped LaGaO3 system as a consequence of the applied biaxial strain. Figure 1 shows that the applied strain produces a net change in unit cell volume, therefore resulting in a nonperfect Poisson behavior for this material. The average Poisson coefficient calculated according to eq 2 is approximately 0.27, in excellent agreement with the experimental values of 0.27−0.29 obtained for this material.28 Here, the Poisson coefficient, ν, is defined by eq 2 ⎛ −2v ⎞ Δa Δc ⎟· =⎜ ⎝ 1 − v ⎠ a0 c0
(2)
where Δc and Δa are the differences between the lattice parameters at a given strain and the corresponding lattice parameter for the unstrained system. A nonperfect Poisson behavior for perovskite-type materials was recently reported for CaMnO3 based on DFT calculations.29 The lattice energy of LaGaO3 as a function of unit cell volume is plotted in Figure 2. Data were fitted to a third-order
Figure 3. Energy-minimized structures for LaGaO3 at −5%, 0%, and +5% applied biaxial strain in the ab plane.
inspection of these figures shows that, as expected, the applied strain in the ab plane has an impact on both the tilting and the mutual rotation of the octahedra. On the basis of the energyminimized structures, bond lengths and angles were calculated and the results plotted as a function of applied strain in the ab plane. Figure 4 shows the evolution with strain of the average Ga− O bond lengths in plane (Ga−Oeq, 4 equatorial Ga−O2 bond lengths) and out of plane (Ga−Oap, 2 apical Ga−O1 bond
Figure 2. Calculated lattice energy as a function of lattice volume for orthorhombic LaGaO3. Fit to a third-order Birch−Murnaghan equation of state (solid line).
Birch−Murnaghan equation of state, and the corresponding calculated bulk modulus is BEPP = 233 GPa. This value is in fair agreement with experimental data obtained for the LaMn1−xGaxO3 system in orthorhombic symmetry, for which the bulk modulus was found to be in the range 153−188 GPa for 0.2 < x < 0.6 and showing an increasing trend along with x.30 The good agreement between calculated and experimental mechanical properties values obtained in this study as well as the reported good agreement between experimental and calculated structural parameters and dielectric constant data
Figure 4. Calculated average bond lengths as a function of strain for LaGaO3. C
dx.doi.org/10.1021/jp509413w | J. Phys. Chem. C XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry C
Article
Figure 5. (a) Calculated average bending angles as a function of strain for LaGaO3. (b) Graphical representation of the bending angles.
3.2. Defect and Transport Properties. Energy-minimized structures at various applied strains were used for calculation of basic defect formation energies, dopant substitution energies, and binding energies for charged and neutral dopant-oxygen vacancy clusters and also to derive migration energies for oxygen transport within the strained and unstrained systems. In the following, the evolution of such quantities with strain will be analyzed in detail. The oxygen vacancies formation energies in the equatorial and apical positions of the GaO6 units are almost equivalent for the unstrained system (Figure 6). These quantities sensibly vary
lengths). Moving from compressive to tensile strain, the GaO6 octahedral environment changes from elongated to compressed along the c axis, with the octahedra belonging to the strain-free model being the more regular with respect to bond lengths. Two kinds of bending angles between linked octahedra were considered for the orthorhombic system, as in our previous work:31 Ga−Oap−Ga (hereafter out of plane), indicating the bending along the c axis, and Ga−Oeq--Ga (hereafter in-plane) as a measure of the bending in the ab plane (Figure 5). Both these angles should be 180° for an ideally cubic perovskite. Bending in the c direction is particularly affected by the imposition of strain in the ab plane and, with a change from 170° to 159° when changing from −5% to +5% strain in the ab plane, is responsible for a major distortion of the structure in this direction when compared to the ideal 180° angle for the perfect cubic perovskite. Compressive strain increases the bending of octahedra in the ab plane whereas sensibly decreasing bending in the c direction with respect to zero strain. The in-plane bending angle presents a maximum at strain values closed to zero, whereas both compressive and tensile strain result in an increased deviation from the ideal value for a nondistorted perovskite structure; the percentage variation of the in-plane bending angle is smaller compared to the out-ofplane bending. Considering the structural data presented in Figures 4 and 5, it can be concluded that the distortion of the unstrained system is mainly due to the bending between adjacent polyhedra, rather than to octahedra elongation along preferential directions. Application of biaxial strain results in an increased distortion of the octahedra due to the differentiation of apical and equatorial bond lengths as well as to a modification of the bending degree between adjacent octahedra. Deviation of the B−O−B angle from 180° in perovskite structures (B being the cation in the octahedral environment) is considered detrimental to the electronic conduction properties of perovskite-type materials, since it lowers the degree of orbital overlapping, so reducing the efficacy of hopping in a small polaron system.32 Even though this is not the case for LaGaO3, which is an electronic insulator, this result may be of more general interest for ABO3 perovskites and relevant, for example, for cobalt-, manganese-, or iron-based perovskite-type oxides used as cathode materials in SOFCs. Indeed, experimental research efforts are now devoted to optimization of materials properties through proper lattice strain in perovskites such as LaCoO3.20,33
Figure 6. Energy difference between the calculated oxygen vacancy formation energy for LaGaO3 at a given strain and the value in condition of zero applied strain for the oxygen equatorial (eq) and apical (ap) positions.
with applied strain, and in particular, for large tensile strains the formation energy of the oxygen vacancy in the equatorial plane is considerably smaller than that in the apical position. As previously found for CaMnO3,29 this result suggests a preference for vacancy ordering in this crystallographic position, which in turn may influence the functional properties of the material. Calculated activation energies for oxide-ion migration in the system are reported in Figure 7a. Two curves are shown: one is related to migration of an oxide-ion vacancy between adjacent positions in the equatorial plane of the GaO6 octahedra, i.e., the D
dx.doi.org/10.1021/jp509413w | J. Phys. Chem. C XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry C
Article
Figure 7. (a) Calculated oxygen vacancy migration energy for LaGaO3 as a function of strain. (b) Graphical representation of the out-of-plane (green) and in-plane (purple) migration paths.
within the range x = y = 0.05 and x = y = 0.20,38 but also values in the range 0.602−0.837 eV were calculated for the same compositional range,39 up to activation energies higher than 1.1 eV for the low-temperature regime (
