Adaptive kinetic Monte Carlo simulation of solid oxide fuel cell components

Ionic conductivities in the solid oxide fuel cell (SOFC) electrolytes yttria-stabilised zirconia (YSZ), calciastabilised zirconia (CSZ), gadolinium-doped ceria (GDC) and samarium-doped ceria (SDC) and the cathode material lanthanum strontium cobalt oxide (LSCO) are directly calculated using DL_AKMC, an adaptive kinetic Monte Carlo (aKMC) program which assumes limited a priori knowledge of the kinetics of systems. The materials were simulated over several milliseconds and over the range of experimentally most relevant temperatures and dopant concentrations (2–18 mol% for doped zirconia, 5–25 mol% for doped ceria and 5–80 mol% for LSCO). Ionic conductivities of the electrolytes at 1000 K are in good agreement with the observed values: CSZ in the range 3 10 3 to 1 10 2 S cm 1 depending on dopant concentration, YSZ 4 10 3 to 3 10 2 S cm , GDC 1 10 2 to 5 10 2 S cm , SDC 1 10 2 to 7 10 2 S cm . LSCO is predicted to have an ionic conductivity of the order of 10 2 to 10 1 S cm 1 depending on Sr content. Average activation energies over all migration processes are 0.4– 0.5 eV for the stabilised zirconias and 0.2–0.3 eV for the doped cerias and 0.3 eV for LSCO, in agreement with experiment. aKMC provides a distinct advantage over traditional KMC methods, in which one has to provide a list of system state transitions. Here, all of the state transitions are dynamically generated, leading to a more accurate simulation of the kinetics as the system evolves.


Introduction
Due to growing concern about global warming and dwindling supplies of fossil fuels, increasing interest is being directed towards alternative sources of energy. Popular and promising alternatives are fuel cells, and a substantial amount of research and development has been undertaken to improve the materials used in these devices. In solid oxide fuel cells (SOFCs) the oxidant (e.g. air, O 2 ) is reduced at the cathode. The oxide ions produced are then transported through the solid electrolyte material, ideally a purely ionic conductor, to the anode where the fuel (e.g. H 2 , hydrocarbons) is oxidised. Electrons from this process then ow from the anode to the cathode, completing the circuit and generating power. Lowering the operating temperature of SOFCs to an intermediate temperature range of 600-800 C is an ongoing area of intense research and many materials have been suggested as suitable electrolytes 1-3 and electrodes [4][5][6][7][8][9] for this purpose.
Stabilised zirconias doped with lower valency ions such as yttrium (YSZ) or calcium (CSZ) along with ceria doped with gadolinium (GDC) and samarium (SDC) exhibit very high ionic conductivities associated with mobile oxygen vacancies formed as a consequence of the doping and are typically used as electrolytes in SOFCs. Previous computational studies of YSZ have tended to use molecular dynamics (MD) 10 or density functional theory (DFT) 11,12 in combination with kinetic Monte Carlo (KMC) to investigate the oxide ion diffusion kinetics, 13,14 as timescales tend to be limited if MD or particularly DFT is used in isolation. Typically common atomistic KMC models employ an on-lattice approximation which limits their ability to describe a system which undergoes large structural changes. This approach also requires a list of possible event mechanisms determined a priori, through experimental and theoretical methods, by estimation or even guessing. This is naturally a severe limitation as the processes involved in the atomic motion are not necessarily intuitive and can be extremely difficult to predict in advance. 15 Furthermore, as the simulation advances the structure of the system changes leading to new possible transitions and altering the activation energies of existing transitions. To overcome these limitations various 'on-the-y' approaches, such as adaptive KMC (aKMC), 16 have been proposed. These are designed for an off-lattice system and the saddle-points between different possible states of the system are located as the simulation progresses. These methods overcome the limitations of the traditional KMC approaches by avoiding the need for a list of mechanisms provided beforehand. State transitions not predicted in advance are allowed, thus permitting a thorough search of the energy landscape and ensuring that the kinetic model is as realistic as possible. As in traditional KMC, certain assumptions are made when using aKMC, namely that the system transitions adhere to harmonic transition state theory. 17,18 Here, we also assume that the preexponential factors, used in the calculation of transition rates, are equal to a typical vibrational frequency (10 13 s À1 ).

aKMC
The key to the aKMC approach is an efficient way to search for saddle points that link the current state of the system to another. One popular approach is the 'dimer' method, 19 the implementation of which used in this work has been described in detail elsewhere. 20 The dimer algorithm begins with a random starting position within the energy basin and progresses by climbing uphill along the lowest eigenvector corresponding to the lowest value eigenvalue of the Hessian matrix, eventually reaching a saddle point at the top. This is a particularly efficient method as only rst derivatives of the energy 21 are required. If all the saddle points leading to different system states can be located, the activation energies of each of these pathways can be supplied to the KMC procedure and, with the pre-Arrhenius frequency factor, the rates of each pathway determined. The system can then be propagated in a dynamically correct way to the next state, and the entire procedure is repeated.
In practice, however, it is computationally demanding (effectively impossible) to nd all the saddle points and demonstrate that all have been found, and the number of saddle points bounding a state grows exponentially with the dimensionality of the system. Xu and Henkelman 22 have developed a method to quantify a 'condence limit' for each state that enough relevant saddle points have been found to progress the simulation with a pre-dened level of condence and this approach has been implemented here. Their method supposes that all of the relevant kinetic events possess activation energies within the range mk B T of the lowest barrier process. m should be carefully chosen, and large enough that relevant transitions are not omitted. For example, for m ¼ 20, an event at the upper limit will be e À20 z 10 À9 as likely to occur as one at the lower limit, assuming that the pre-factors are equal. Events that occur within these limits are therefore relevant to the kinetics of the system, and the condence (C) that a relevant saddle has not been missed is: N r is the number of sequential searches that nd relevant, but redundant (non-unique) processes and a lies between zero and one and describes the relative probability of nding each relevant saddle point. a ¼ 1 describes a system where there is an equal probability of nding each relevant saddle point, and is the value used in this work. More justication of eqn (1) is detailed elsewhere. 22 By setting a condence limit of, e.g. 95%, the simulation will run until N r ¼ 20 searches nish without nding a new, unique saddle point. Each saddle point search proceeds assuming no knowledge of the local environment, with no saddle point recycling.
We use a new adaptive KMC program, DL_AKMC, † as follows. Initially the provided material structure is minimised using a user-specied algorithm as implemented in DL_FIND. 23 Possibilities are steepest descent, conjugate gradient following Polak-Ribiére, 24 L-BFGS, 25,26 P-RFO, [27][28][29][30] Newton-Raphson/ quasi-Newton, damped dynamics, random (stochastic) search 31,32 or by a genetic algorithm. [33][34][35] This minimised structure is then used as the initial basin for the task-farmed saddlepoint searches either by the improved dimer method 19,20 or the NEB method 36 as implemented in DL_FIND, as specied by the user. These dimer searches are initiated near the local minimum by displacing the system away from the minimum in a random direction. The algorithm for displacement used here is generaleach atom is displaced up to a maximum of AE0.4Å.
The kinetic Monte Carlo algorithm has been documented elsewhere 37,38 and Fig. 1 gives a brief overview of the parallelisation used. Once all of the transitions have been identied, the rate n i for each event i is then obtained from eqn (2): n 0 is the pre-exponential factor, here set equal to a typical vibrational frequency (10 13 s À1 ) as is usual in the absence of more detailed information. 12 DE is the activation energy barrier for ion migration determined using the dimer method, k B the Boltzmann constant and T the temperature. One transition is chosen at random with a probability proportional to the relative rate of the transition. The timestep of the simulation is then advanced by Dt (eqn (3)): Fig. 1 Schematic of the parallelisation and simulation process implemented in DL_AKMC.
u is a random number between zero and one. The system structure is then updated to reect the transition chosen, and the saddle-point searching algorithm begins again. This process continues until the maximum simulation or computational time is reached.

Ionic conductivity
Key to the operation of an SOFC electrolyte or cathode is its ionic conductivity. For efficient functioning of the fuel cell both electrolytes and cathodes should possess high such conductivity. This is related to the oxygen self diffusion coefficient by the Nernst-Einstein relation (eqn (4)): s i is the ionic contribution to the conductivity from species i, z i e the charge of species i, c i the concentration of ionic defects (in this paper, oxide-ion vacancies) and D i the diffusion coefficient of species i. This approach is strictly only valid for dilute systems, and not our doped systems which are almost always concentrated solid solutions. Grope et al. 39 and Pornprasertsuk et al. 11 have simulated such systems circumventing this issue by direct inclusion of an applied eld. In their approach, the electrolyte is split into a series of 'slices' that are innite in two dimensions (y-and z-) perpendicular to the direction of the applied eld (see Fig. 2). The applied eld alters the activation energy barrier for ion migration (eqn (5)): DE corr is the new, corrected activation energy barrier for ion migration, DE 0 the activation energy barrier obtained from the saddle-point searches, V shi is the total potential difference between the destination slice and the initial slice in which the vacancy is located. a is termed the 'symmetry factor' and has a value between zero and one which is dependent on the prole of the migration barrier under the applied potential but in this work is set to 0.5. For any oxide-ion vacancy movement that occurs within a slice, V shi is zero. The total potential, V i total , at each slice, i, is split into two termsthe electrode potential (V i elec ) and the space charge potential (V i sc ): V i elec is the potential on slice i due to charge accumulation at the electrodes. Electrons are assumed to have a much higher mobility than oxide ions and so V i elec is equal to the difference between V i total and V i sc . In our simulations, V elec is assumed to vary linearly across the electrolyte: V 0 elec , the electrode potential in the initial slice, labelled 0, is set to zero for convenience. V N elec is the electrode potential on slice N, the nal slice in our conguration. The space charge potential per slice, i is given by: V 0 sc , the space charge potential in the initial slice, 0, set to zero. The slice spacing of all the cubic unit cells of the crystal structures considered in this work is half the lattice parameter, a. In using eqn (8) during the simulation, ionic movement is assumed not to change the geometry of the simulation cell. E i is the electric eld associated with each innite slice, i given by: 3 0 is the vacuum permittivity and 3 r the relative permittivity of the material, set at 40 in our simulations, equal to the experimental value for 8 mol% YSZ. 40 Small variations in the permittivity that occur when simulating different materials do not greatly inuence the results. The charge density per slice (r i ) is the total defect charge in the slice divided by the surface area of the periodic slice. To simulate an impedance measurement an alternating potential, V appl , is applied: t is the time, f the frequency and V 0 the amplitude of the applied potential. V N total (eqn (6)) is set to ÀV appl . Under the applied eld the activation energies and subsequently the rates (DE and n i respectively in eqn (2)) for migration of charged species in the eld direction change, and the resulting displacement can be used to calculate the contribution to conductivity from species i, s i (eqn (11)): hXi is the mean displacement of species i in the applied eld, E eld , c X the concentration of i and t is the simulation time. We assume no change in charge as the charged species migrate.
For the aKMC calculations the materials were simulated for a minimum of 1 ms, and a summary of the parameters required is given in Table 1.

Application to oxides
We have investigated four solid oxide fuel cell (SOFC) electrolyte materials with the uorite structure: yttria-and calcia-stabilised zirconica (YSZ and CSZ respectively, with dopant concentrations varying from 2-18 mol%) and gadolinium-and samarium doped ceria (GDC and SDC respectively, with dopant concentrations varying from 5-25 mol%), and one cathode material with the perovskite structure, strontium-doped LaCoO 3 (LSCO, with dopant concentrations varying from 5-80 mol%). Pictures of the two structures are shown in Fig. 3. YSZ, CSZ, SDC and GDC are all stable at intermediate temperature SOFC operating temperatures (up to 900 C), and YSZ, SDC and GDC all possesses high ionic conductivities ($10 À2 S cm À1 at 700 C). 41,42 Doping of zirconia involves substitution of some some Zr 4+ ions with Y 3+ or Ca 2+ and of ceria replacement of some Ce 4+ with Gd 3+ or Sm 3+ . To maintain charge neutrality, one oxygen vacancy is created for every two Zr or Ce ions replaced with Y, Gd or Sm ions. When Zr ions are replaced with Ca ions, one oxygen vacancy is created per ion replaced. This gives rise to a considerable number of oxide ion vacancies which greatly increase the ionic conductivity of these materials relative to the undoped systems.

Potential model
For the adaptive kinetic Monte Carlo simulations a rigid-ion model, including 2-body Buckingham potentials is used (eqn (12)) to model the atomic interactions. The potential parameter set is listed in Table 2. A non-Coulombic potential cutoff of 12Å was used throughout. Each ion is assigned an integer charge, i.e. +2 for Ca and Sr, +3 for Gd, Sm, La and Co, +4 for Zr and Ce and À2 for O.
A, r and C are constants, r ij is the distance between ions of type i and j, and z i e and z j e are the charges of species i and j respectively. Since our potential set uses a number of potentials from different sources, we have validated the set by checking good agreement between simulated and experimental values for bulk properties (lattice parameter, elastic constants and bulk modulus) of a number of binary oxides using the GULP code. 50 The results are shown in the ESI, Section 3. † Cubic supercells containing 4 Â 4 Â 4 cubic unit cells (for SDC, GDC, YSZ and CSZ) and 6 Â 6 Â 6 cubic unit cells (for LSCO) were constructed and both the dopant cations and required oxygen vacancies were placed randomly in the lattice, subject to the constraint that each slice (0 to N, see Fig. 2) has zero net charge. It has been suggested that the oxygen vacancies in YSZ preferentially occupy sites near to yttrium ions (second nearest neighbour positions 51 ), which is not reected in our choice of initial structure (randomly distributed stabilisers and vacancies). We have performed studies investigating the inuence of the starting structure on the ionic conductivities at the end of the adaptive kinetic Monte Carlo simulations. Both structures with anion vacancies and dopant cations as nearest neighbours, and structures with anion vacancies and dopant cations separated by a minimum of 4.5Å result in ionic   conductivities similar to those seen with the random arrangement used here. The initial arrangement of the dopant atoms and vacancies appears to bear little inuence on the ionic conductivity of the material at the end of the simulations. More detail on this is given in the ESI, Section 1. † The initial random structures used in this work were equilibrated for 20 ps (NPT ensemble) using DL_POLY 52 at the temperature used in the subsequent DL_AKMC simulation (typically 300 K, with the exception of those used to create the Arrhenius plot in Fig. 9 later). Aer equilibration each slice was checked to ensure charge neutrality was maintained, with the resulting arrangements used as input structures for DL_AKMC.
To check further the quality of the potentials employed, the lattice parameters of the simulation cells (aer equilibration) are compared with experimental lattice parameters as a function of the dopant concentration ( Fig. 4 and 5). Fig. 4 shows that the calculated lattice parameters for the doped zirconias and cerias agree with the experimental values with a match of <1%. Similarly, the lattice parameters for LSCO are in good agreement, with the calculated lattice parameter slightly larger then experiment (by 0.6-1.6%). Both experiment and simulation agree as to the increase of the lattice parameter of the doped system with increasing dopant content. Our potentials thus describe these systems adequately and are suitable for use with adaptive kinetic Monte Carlo.

Results and discussion
Conductivities were calculated (eqn (11)) for all the materials and the results are shown in Fig. 6 and 7 for the electrolyte materials and Fig. 8 for LSCO, along with some representative experimental values. It can be seen that there is good agreement between the simulations and experiment for all of the materials modelled. In particular the calculated values of the ionic conductivities in the doped cerias (Fig. 7) match very closely the experimental magnitude of the conductivity, and its variation with dopant concentration. The calculated maximum in conductivity for SDC occurs at 16 mol%, higher than that found experimentally by 2 mol%. Our choice of electrolytes has been governed by their high ionic conductivity and low electronic conductivity and so we would expect that our simulations, which solely determine the ionic conductivity, should agree with the experimental values which are the sum of ionic and electronic contributions. The good agreement found between the calculated and experimental conductivities is emphasised by our use of linear y-axes, in contrast to the logarithmic scale typically deployed for such gures.
The form of the ionic conductivity versus dopant concentration gures has been discussed by others. 55,56 For YSZ in particular, the decrease in conductivity at larger dopant concentrations has been attributed, at least in part, to the formation of strongly-bound Y 0 Zr complexes reducing the availability of free V $$ O available for diffusion. 57 CSZ exhibits a more unusual curve, with two maxima appearing in   This journal is © The Royal Society of Chemistry 2014 the dopant range simulated. While the location and magnitude of the second maximum at 14 mol% matches well with the experimental data available, 58 the general paucity of experimental ionic conductivity measurements for CSZ is such that there are no data for comparison at 6 mol%, and so the rst maximum remains an interesting prediction for future verication. It is possible this rst maximum may be simply an artifact of the simulation, and a discussion on the reproducibility of results and inuence of the starting conguration is given in the ESI, Sections 1 and 2. † For a mixed electronic and ionic conductor such as LSCO the electronic conductivity s el $ 10 2 -10 3 S cm À1 , dominates. 59 The simulated values of the ionic conductivity (s ion ) shown in Fig. 8 cannot be compared directly with experiment. However values of s ion have been calculated from experimental oxygen diffusion coefficients and the Nernst-Einstein equation. These values are in the range 10 À3 to 10 À1 S cm À1 at 900 C depending on composition. 59 These experimentally-derived values, relate to different temperatures (1048-1173 K), but do highlight the broad margin of error inherent in these measurements. Our theoretically determined values fall within the experimental range of ionic conductivities, and exhibit a coherent trend across the entire dopant range. A maximum in the ionic conductivity is found at $48 mol% Sr. Combined with the good agreement between theory and experiment for the electrolytes, we are condent that our values for the ionic conductivity of the LSCO cathode are accurate. Materials such as LSCO, where experimental values are difficult to obtain, emphasise the importance of these simulations for extracting important physical properties that would otherwise be difficult or impossible to measure experimentally.
For each material, one experimentally common dopant concentration was chosen (10 mol% for YSZ and CSZ, 20 mol% for SDC and GDC and 50 mol% for LSCO) and simulated over a temperature range of 600-1100 K. The conductivities were determined over these temperature ranges and the resultant Arrhenius plot is shown in Fig. 9. The average activation energy of the migration processes occurring in each material can be determined from the gradient of a linear t to the data.     9 Arrhenius plots for ionic conductivities of 10 mol% YSZ (black, circles), 10 mol% CSZ (red, squares), 20 mol% GDC (green, diamonds), 20 mol% SDC (blue, triangles) and 50 mol% LSCO (orange, crosses). content, with the activation energy dropping to 1.4 eV for 50 mol% Sr from 2.2 eV for 20 mol%. In general we nd that the calculated activation energies for the migration processes in these materials are lower than those found experimentally. This can possibly be attributed to the use of a xed charge, rigid-ion model where the ionic charge does not change throughout the simulation; in reality the charges on the ions are likely to change depending on the local environment.

Conclusion
The reliability of the adaptive kinetic Monte Carlo program DL_AKMC has been veried by the accurate simulation of ionic conductivities of the solid electrolytes yttria-and calcia-stabilised zirconica (YSZ and CSZ respectively), gadolinium-and samarium doped ceria (GDC and SDC respectively), and of one cathode material with the perovskite structure, strontiumdoped LaCoO 3 (LSCO). The calculated ionic conductivities of all the materials are of the same order of magnitude as those found experimentally. The variation of ionic conductivity with dopant concentration follows the same form as experiment. The simulations also indicate the presence of a possible peak in the conductivity of CSZ at $6 mol%, which would be of interest to probe experimentally. Average activation energies of all migration processes are between 0.4-0.5 eV for the stabilised zirconias, between 0.2 and 0.3 eV for the doped cerias and 0.3 eV for LSCO.
This work has showcased the ability of the program to simulate experimentally measurable and important properties, while retaining the exibility of using common potential types and the ability to determine transition states on-the-y in a massively parallel fashion. Measuring solely the ionic conductivity in mixed ionic-electronic conductors is difficult, and there is a large variation in the experimental ionic conductivity, particularly with LSCO. The ability for a simulation to model reliably the ionic conductivity of such mixed conductors with high precision and accuracy is of critical importance, and such predictions are now achievable and efficient. Given the exibility and general applicability of the DL_AKMC program, it will likely prove to be a vital tool in a wide range of scientic elds.