Modeling of complex interfaces: Gadolinium-doped ceria in contact with yttria-stabilized zirconia

Gadolinium-doped ceria (GDC) and yttria-stabilized zirconia (YSZ) are well-known electrolyte materials in solid oxide fuel cells (SOFCs). Although they can be used independently, it is common to find them in combination in SOFCs, where they are used as protective layers against the formation of secondary phases or electron conduction blockers. Despite their different optimum operating temper-atures, it appears that oxygen conduction is not affected by their interface. How-ever, the intrinsic mechanisms of oxygen diffusion at these interfaces still remain unclear. One of the main difficulties when modeling the contact between different materials, or indeed different particles of the same material, is caused by the structural complexity of these systems. If we wish to evaluate the properties of the materials, we first need to obtain a model that includes the main features of the GDC/YSZ interface, such as large-scale defects or cation interdiffusion in the contiguous phase. Since the generation of such a mixed system is complicated, we show here how the “ amorphization and recrystallization ” strategy can help us to obtain realistic systems. In this, the first of our papers on the structure and properties of layered GDC/YSZ materials, we discuss the structural features of the grain boundary between GDC and YSZ obtained by molecular dynamics simulations.


| INTRODUCTION
Solid oxide fuel cells (SOFC) are of significant interest to the scientific community due to their high energy conversion efficiency and fuel flexibility, in combination with a very low environmental impact. 1 SOFCs comprise three main parts: the anode, which normally consists of metal nanoparticles supported on a metal oxide that also acts as the electrolyte; the electrolyte itself; and the cathode, which is normally a perovskite material.
At the anode, CO and H 2 react with the O 2À that is generated at the cathode from the reduction of O 2 and transported through the electrolyte, producing CO 2 and H 2 O. Normally, anodic reactions take place at the triple phase boundary (TPB), which is where the reaction mixture, the metal nanoparticles, and the electrolyte meet. 2 Thus, after the reaction takes place, the oxygen can be automatically driven toward the cathode. Gadolinium-doped ceria (GDC) and yttria-stabilized zirconia (YSZ) are two wellknown doped metal oxides used for that purpose. [3][4][5] Although isostructural, they exhibit different properties, and one of the more important differences is in their working temperatures. The required temperature for optimal ionic conduction in YSZ is between 873 and 1273 K, whereas in GDC this temperature is lower, that is, between 775 and 1075 K. 6,7 This variance in optimum working temperature can be explained in terms of the different contributions to the activation energy of the oxygen diffusion process, which are the association and the migration energy. The former relates to the energetic stabilization of oxygen vacancies in the materials due to their Coulombic interaction with the gadolinium and yttrium dopants in GDC and YSZ, respectively, whereas the latter describes the energy required to move one oxygen to a neighboring vacancy. Association energies for ½Gd 0 Ce À V ÁÁ O À Gd 0 Ce are slightly lower than for ½Y 0 Zr À V ÁÁ O À Y 0 Zr , that is, around 0.10 and 0.20 eV for GDC, but between 0.3 and 0.35 eV for YSZ. However, the differences between the migration energies are larger, with GDC showing values between 0.70 and 0.80 eV while YSZ has migration energies between 0.90 and 1.25 eV, thus showing that global activation energies for GDC are easier to overcome than for YSZ. 8,9 Unfortunately, despite their extensive usage, both materials present some drawbacks. YSZ has been shown to react with the cathode material, which is normally La (1Àx) Sr x MnO (3Àx/2) (LSMO) or similar perovskites. As a consequence, highly resistive phases of SrZrO 3 or La 2 Zr 2 O 7 are formed, reducing the performance of the cell. 10,11 One way to prevent their formation is to use a protective layer between the materials, which is normally GDC. [12][13][14] Hence, it would seem logical to use GDC instead of YSZ, but unfortunately above 600°C and under low oxygen partial pressure, Ce 4+ is reduced to Ce 3+ , thus inducing electronic conductivity. The electronic conductivity can be blocked while maintaining the ionic conductivity through the introduction of a thin layer of YSZ in the GDC matrix, 15 this GDC/YSZ bilayered electrolyte also presents some inconveniences. During sintering, GDC and YSZ can react, potentially leading to the formation of an interdiffusion layer, which can also decrease the performance of the fuel cell. 16,17 Although mixtures of the two materials are used as electrolytes, little is known about the interface between GDC and YSZ. Even though in-depth understanding of the reactivity and its structural consequences on the materials, as well as the diffusion behavior of the oxygen anions, may suggest improvements to the materials' performance. There are several techniques that can be used to model the interface, but the easiest and most practical one is by molecular dynamics simulations, which has been shown to produce accurate and reliable results for related materials. [18][19][20] Here, we therefore introduce a strategy that is capable of modeling, in a realistic fashion, the interface between GDC and YSZ, using interatomic potentials and molecular dynamics simulations.

| COMPUTATIONAL DETAILS
The generation of a realistic model for the interface between gadolinium-doped ceria (GDC) and yttria-stabilized zirconia (YSZ) was achieved using a combination of the METADISE code 21,22 and the DL_POLYv4 code. [23][24][25] Both codes used the same types of interatomic potentials to describe the interaction between the different atomic species, with the main difference that the META-DISE is used only for static calculations, whereas DL_POLY is used for the molecular dynamics simulations. A brief description of the theoretical background to METADISE, and how it treats surfaces and interfaces can be found in the electronic supporting information (ESI) and Figure S1.
Since both materials are ionic solids, we have assumed the Born model in which the interaction between ions is described by attractive long-range Coulombic terms, balanced by short-range repulsive interactions and van der Waals forces. These latter interactions are evaluated in both codes using the Buckingham potential, expressed as in Equation 1: where A ij (eV), q ij ( A), and C ij (eVÁ A 6 ) are parameters between atoms i and j that are used to fit the equation to experimental data. The Buckingham potentials used throughout this work are listed in Table 1, [26][27][28] and they all have proved to reproduce lattice parameters, volume expansion, or defect association energies for the corresponding materials GDC and YSZ. [29][30][31][32][33] The electrostatic interactions are evaluated in both programs using the Ewald summation, which calculates this interaction in both the real and the reciprocal space, also including a self-interaction correction term. It is worth noting that DL_POLY4 uses the smoothed particle mesh Ewald (SPME) method, 34 which, however, only differs from the standard Ewald summation in the treatment of the reciprocal space terms. Finally, two different ensembles were used in the molecular dynamics simulations: the microcanonical ensemble (NVE) and the Hoover isobaric ensemble (NPT). 35 T A B L E 1 Buckingham potential parameters for the different anion-anion and anion-cation interactions. Cation-cation interactions are assumed to be 0 The F between these two materials can thus be calculated as: where a CeO 2 and a c-ZrO 2 are the lattice parameters of each material, and n and m are the number of repetitions of each unit cell, respectively. In Equation 2, n relates to CeO 2 and m relates to c-ZrO 2 . Hence, using the aforementioned bulk parameters and Equation 2, F between CeO 2 and c-ZrO 2 is 6.01%. According to the experimental lattice parameters for YSZ and GDC, which are 5.423 A and 5.140 A, respectively, 39,40 the lattice misfit that we would expect for them is 5.36%, and it has been reported to be around 6.00% for the GDC and YSZ bulk materials, depending on the dopant concentration. From experimental data, we know that this lattice misfit, at least in the case of CeO 2 /YSZ, is accommodated through dislocations in both materials, 41,42 so we would expect similar structural defects for the GDC/YSZ interface. Another factor to bear in mind is cation diffusion. Harris et al. 43 indicated that, when considering heteroepitaxial systems, cation exchange at the interface is likely to occur, generating a so-called interdiffusion layer. In view of these structural complexities, the construction "by hand" of a realistic model becomes impossible.
The most plausible strategy to obtain a model that includes all the aforementioned defects is the amorphization and recrystallization (A&R) strategy, which has been pioneered and used by Sayle and coworkers in a number of different systems, including SrO/MgO 20 and CeO 2 / YSZ. 31,33 This strategy can be divided into two different steps. First, we generate an ideal structure that minimizes the lattice misfit between GDC and YSZ, by finding a specific near coincidence site lattice (NSCL). 33 Next, we search for a lower energy configuration, using the previous system as a starting point, by increasing the ionic mobility, that is, by bringing the system to a temperature close to the melting point. In time, we recrystallize the system by applying a high pressure while maintaining a high temperature. Finally, once the system is recrystallized, we release the pressure and gradually cool it down until we reach a temperature close to 0 K.

| Near coincidence site lattice
The near coincidence site lattice (NCSL) consists of finding particular supercells for both GDC and YSZ, which minimize the lattice misfit between the two materials as much as possible. Normally, the lattice parameters that are used are those belonging to the respective bulks. However, since the (111) surface is the facet that is the contact between GDC and YSZ, 44,45 we decided to use their surface lattice parameters instead. Accordingly, it seems logic to use the experimental lattice parameters for GDC and YSZ to find the corresponding NCSL. However, we decided to use the calculated lattice parameters for CeO 2 and ZrO 2 instead, mainly because of two reasons. First, any model is usually based on the calculated lattice parameters and not on the experimental ones. Consequently, we should calculate the GDC and YSZ parameters, but that would involve calculating a large amount of unit cells with different dopant concentrations and distributions. Second, if we were using the doped bulks, the surfaces and grain boundaries that we would obtain would suffer a great reconstruction, because of the impossibility of finding nonpolar terminations.
We thus started by generating a slab model for the (111) surface of both parent materials (CeO 2 and c-ZrO 2 ) using the METADISE code. The (111) surface for a fluorite-type system has two different nonpolar terminations: metal (M t ) and oxygen (O t ), as shown in Figure 1B,C, respectively. The relaxed surface energy (c r ) obtained for these terminations was 3.90 and 1.58 JÁm À2 , respectively, indicating that the O t is the most stable termination, as expected from the fact that it is a nonreconstructed stoichiometric Tasker type II surface, 22 whereas M t is a reconstructed dipolar type III surface, which are generally unstable surfaces.
Next, the two surfaces were used to generate the grain boundaries (GBs). There are two main types of GBs: tilt and twist. Tilt GBs are obtained when the rotation axis of one grain is parallel to the GB's plane, whereas twist GBs are obtained when this rotation axis is perpendicular to the GB's plane (see Figure S2 in the ESI). In this work we only considered tilt GBs generated by reflection of the slab model. Since c r relates to the stability of the surface but not to the GB, we studied the GBs originated by reflection of both M t and O t terminations and calculated the relative stabilities of the two GB systems. As c-ZrO 2 is not stable under standard conditions, 46 we performed this step only for CeO 2 , assuming that the results we obtained would be transferable to c-ZrO 2 .
After slab reflection, we calculated the interface potential energy surface (IPES) by scanning the relative position of one grain with respect to the other. This procedure allowed us to identify the grain arrangement that lowered the energy for both terminations, as shown in Figure 2. For simplicity, from now on we will label these grain boundaries GB-O t and GB-M t , respectively.
For the lowest energy GB arrangement, we calculated the formation energy (E f ) and the cleavage energy (E c ), according to Equations 3 and 4: where E GB is the total energy of the grain boundary, E B is the bulk energy, and E S is the energy of the surfaces. According to Equations 3 and 4, E f is the formation energy with respect to the bulk materials and E c accounts for the energy required to separate two grains thereby exposing their surfaces. Therefore, the most favorable GB will be the one that minimizes E f but maximizes E c . E f for both GBs is calculated to be approximately the same: 1.87 JÁm À2 for GB-O t and 1.94 JÁm À2 for GB-M t , respectively. Interestingly, E c for GB-O t is 1.29 JÁm À2 , but for GB-M t is 6.26 JÁm À2 , that is, despite having approximately the same formation energy, GB-M t is more stable than GB-O t , conversely to what we observed with the clean surfaces. Consequently, we have continued our study focusing only on the M t grain boundary.
In the M t slab, we next substituted certain Ce and Zr for Gd and Y and generated the corresponding oxygen vacancies in order to keep the electroneutrality of the system. The dopant concentration is derived from the general stoichiometric expression of GDC, which is Ce (1Àx) Gd x O (2Àx/2) , which also applies for YSZ as . Hence, our model simulates a dopant concentration of x=0.14, that is, 14% for both GDC and YSZ. Although this concentration is somewhat higher than the usual experimental one, which are 8% for YSZ and 20% for GDC, it is still found within the range of the acceptable concentrations. 11,16,[47][48][49] In addition, we transformed the hexagonal (111) surface unit cell into a rectangular (111) surface unit cell for computational simplicity. Therefore, we had to locate the NSCL for b and c vectors, respectively, which are both lateral to the surfaces.
Thus, according to Equation 2, n and m are equal to 1, F was 6.78% in both directions, which was minimized to 0.18% for orthorhombic b direction and to 0.17% for the orthorhombic c direction, respectively, by considering n=14 and m=15 in both directions (see Table 2). Hence, the dimensions of the simulation box that will contain both thin layers will be b % 107.70 A and c % 186.50 A. Finally, normal to the surface ã 50.00 A, resulting from the height of the GDC slab of 22.959 A, the YSZ slab of 20.539 A, and an initial distance between the materials, which we set at 3 A per side. If we had assumed the experimental lattice parameters for GDC and YSZ instead, which are listed in Table 2 as well, the lattice misfit would have also been minimized by considering n=14 and m=15, but with a slightly larger F of 1.54%. This difference between lattice misfits is indicative that it is likely that we are potentially introducing strain in our system. If that were the case, we would expect the presence of large-scale defects once the A&R is finished.

| Amorphization and recrystallization
For the A&R, we have used two different ensembles: NVE and NPT. As it has been shown that with constant pressure ensembles, the system retains its original structure, 33 NVE was used for the amorphization step, which consisted of melting both materials for a very short period of time  Table 3. After randomization, we applied 200 katm using an NPT ensemble for to recrystallize the system, which also involved a compression of the cell. The selection of this pressure is independent of any experimental pressure and its only purpose is to speed up the crystallization phase. 33 After 3.3 ns, the system was fully crystallized and we proceeded first to release the pressure, with a production run of 3 ns, followed by a controlled cooling down of the system to 10 K.

| Structural analysis
After applying the A&R process, the dimensions of the system have changed significantly. The height of the unit    From a visual analysis of Figure 3, the bulk does not show any void regions, and full crystallinity seems to be recovered. However, the layer structure is clearly not as regular as it was prior to A&R. In the initial geometry ( Figure 3A), both GDC and YSZ had seven metal atomic layers, but after A&R ( Figure 3B) each material has now between 10 and 12 metal layers, depending on the region. This is a clear indication that the bulk is now less regular than it was before, so the presence of large-scale defects, such as dislocations, is expected.
For a more detailed analysis, we have divided the system into different sections in order to identify the presence of any defects. As a representative sample, we will discuss two different slices, shown in Figure 4. Figure 4A illustrates a section of the material where full crystallinity is almost recovered. However, some local disorder can be observed at the interface between GDC and YSZ, as well as a degree of cation diffusion that has taken place between the two phases. It is also worth noting that the upper GDC phase is slightly bent, which is an expected distortion when there is a small lattice misfit between materials, 41 although no clear dislocations are observed.
In contrast, the section depicted in Figure 4B shows a different scenario. Here, despite the fact that a crystalline arrangement can also be identified, the upper interface is clearly irregular, which agrees with experimental evidence for a roughness between materials. 31 Moreover, we observe a less dense region in the left-hand side of the YSZ phase along the ã direction. Conversely, the GDC shows a major densification in the same region. A closer look reveals that this higher densification is actually a series of successive screw dislocations that are depicted in more detail in Figure 4C. This series of dislocations, primarily located in the upper GDC phase, lead to what appears to be the start of a crack in the lower GDC phase, which is being filled by YSZ, and could explain why the YSZ is less dense in the upper grain boundary. The presence of the dislocations and the slight bending of the GDC phase give rise to apparent different domains, highlighted in blue on the left in Figure 5. It thus seems quite clear that despite the low lattice misfit achieved by using the NSCL based on the CeO 2 and ZrO 2 lattice parameters, this apparent epitaxy was lost during the A&R, and is now accommodated through a series of defects. Also, based on the different distribution of defects across our model, it seems clear that the strain is not uniform within our material, as it happens in any real system.
The presence of these large-scale defects in our model somehow differs from a similar interface obtained using the A&R methodology. Sayle et al. use the A&R to obtain a realistic grain boundary for CeO 2 and YSZ. 33 Unlike us, during the recrystallization phase they were able to dissolve the different grains that were formed during the process. Although we have used the exact same strategy, the main difference between their model and ours is that they base the construction of the NCSL on the experimental lattice parameters of CeO 2 and YSZ.
That unequivocally brings us back to whether our assumption of considering the lattice parameters of the calculated parent materials (CeO 2 and ZrO 2 ) was right or not. As we were expecting, the small strain originated from the remaining s lattice misfit in our NCSL is the responsible of the different large-scale defects that we obtain. Therefore, that implies that unless we are considering the experimental (or calculated) lattice parameters for the doped systems, we will have large-scale defects in our model. However, experimental interfaces do show large-scale defects, which are included in our model. Hence, we conclude that despite not fully minimizing the lattice misfit with our approximation, the overall result is improved by their presence making our system to be more realistic.

| Dopant distribution and clustering
The cation diffusion shown in Figure 4 is corroborated by the atomic percentage profile, illustrated on the right in Figure 5. This profile clearly shows cation penetration between GDC and YSZ, and also reveals that this effect is larger for Ce 4+ and Zr 4+ than for the dopants Gd 3+ and Y 3+ . According to the concentration profile, Ce 4+ and Zr 4+ penetrate about 7 A inside the contiguous phase, and so does Y 3+ , which is the same behavior as observed  50 In contrast, Gd 3+ barely penetrates inside YSZ phase, which is probably a consequence of its initial distribution in the CeO 2 matrix. One could argue that the gadoliniums' atomic mass is related to its lower mobility during the amorphization phase and thus explain why Gd 3+ shows less penetration. However, if true, we would expect a correlation between penetration and atomic mass, and this is not the case.
Cation penetration in the contiguous phase is an expected phenomenon, as the formation of Ce a Gd b Y c Zr e O 2-f mixtures is very common, and among other applications, is exploited to prevent delamination between GDC/ YSZ. 16,51,52 In this case, however, no interdiffusion layer is formed, since the time spent to randomize the system during A&R was deliberately kept short in order to avoid this behavior. Nonetheless, it is worth mentioning that the cation interdiffusion is mainly controlled by the length of the amorphization phase. Therefore, if we were intending to study the formation and behavior of the interdiffusion layer, we would only need to run a longer randomization simulation. Finally, in the concentration profile, we also included the oxygen concentration, which remains constant throughout the different materials.
Aside from the cation penetration, we also assessed whether dopants are uniformly distributed throughout the bulk, which is important in order to evaluate the possible formation of different phases that could develop due to dopant clustering or the formation of an interdiffusion layer. We therefore analyzed the radial distribution function (RDF), as shown in Figure 6.
RDF provides information about the average interatomic distances. So far, since we do not have experimental distances for the doped materials, we are going to use the pure bulk ones. The experimental Ce-O distance is around 2.34 A (based on the 5.411 A experimental lattice parameter), but the RDF indicates that the average Ce-O distance in our system is 2.28 A, suggesting that the GDC phase could be under compression. This finding is consistent with the presence of screw dislocations and the slight bending of the GDC phase, as already discussed. Interestingly, according to the RDF, the Gd-O distances are 2.31 A. The Gd-O distance in the experimental Gd 2 O 3 unit cell oscillates between 2.29 A and 2.32 A, which fit perfectly the RDF of our mixed system. These values, however, do not indicate that the strain is actually influencing the Gd-O distances, but suggest that vacancies and dopants are associated. The presence of an oxygen vacancy next to the gadolinium allows it to be better accommodated in the strained matrix and consequently, Gd-O distances are more similar to the experimental ones. Conversely, this effect could not be seen for Ce-O, because the proportion of Ce atoms and oxygen vacancies is much smaller than for Gd/ V O . It is worth remembering that this structure is obtained at the end of A&R, when the simulation is performed at 10 K, thus the association between dopants and oxygen vacancies is expected.
Regarding the YSZ phase, the Zr 4+ -O 2À distance derived from the experimental bulk parameter is 2.20 A, 37 which is 0.17 A larger than the one we obtained through the RDF analysis, which on average is only 2.03 A. This discrepancy could be an effect of the major irregularities found in the YSZ, with more amorphous regions than in the GDC phase. For Y 3+ -O 2À , experimental Y 2 O 3 shows an average distance of 2.29 A, 53 and our RDF is in excellent agreement with a distance of 2.30 A.
These cation-oxygen distances, however, do not provide any indication whether the dopants are evenly distributed or not. However, the RDF peaks further out from the first neighboring sphere provide this information. The superposition of the peaks relating to the Ce 4+ -O 2À and Gd 3+ -O 2À pair, and the Zr 4+ -O 2À and Y 3+ -O 2À show that the dopants are well dispersed in their respective matrices (ceria and zirconia) and do not form Gd 2 O 3 or Y 2 O 3 domains. 33

| Induced strain
Although dislocations help to accommodate lattice misfit, it is known that at the interface the system is strained trying to accommodate the different lattice parameters. 54 We have shown how the strain can affect the materials in Figure 4A, where the GDC phase is slightly bent. Thus, materials can also become compressed or expanded at the interface, which is one of the factors that potentially leads to fractures. In order to gain further insight into this effect, we have performed A&R calculations, assuming that the entire system was either GDC or YSZ, and compared these pure systems with the GDC-YSZ interfacial system. In this way, we could evaluate the degree of compression or expansion that the materials suffered in the mixed system. From the unit cell values obtained for each material at 10 K after applying A&R, we have calculated the volume variation (DV%) according to Equation 5, where V GDC/YSZ represents the volume of the GDC/YSZ unit cell, and V i represents the volume of either the GDC or the YSZ unit cell. Hence, positive values imply expansion and negative values compression of the interface with respect to the isolated materials.
As observed in the data collected in Table 4, both lattice vectors and volume differences are different depending on the system that we are considering. GDC, when it is in contact with YSZ, suffers a volume compression of 7.91%, whereas YSZ experiences a slightly larger expansion of 8.85%. The lattice vectors behave in the same way, where we observe isotropic compression for GDC and isotropic expansion for YSZ. This compression/expansion pattern, in addition to the screw dislocations and some void regions near the grain boundaries between GDC/YSZ, could potentially have a strong impact on the oxygen conductivity of the system. It is worth mentioning that this effect is intrinsic only to the region near the boundary between GDC/YSZ, and it dissipates as soon as we move away from the interface. 41

| CONCLUSION
In this paper we have used molecular dynamics simulations to investigate the interface between two well-known doped metal oxide materials used as electrolytes in SOFCs: GDC and YSZ. The degree of complexity these interfaces show experimentally requires a sophisticated model that includes all the different types of defects in order to obtain realistic results. To do so, we have used the A&R strategy, 20 using an initial interface that minimizes the lattice misfit between materials according to the near coincidence site lattice protocol.
The resulting mixed material had lost its initial structure and undergone a strong layer reorganization during the crystallization step. As a consequence, the lattice misfit was accommodated via dislocations that were spread through both GDC and YSZ, reflecting a more realistic system than the initial perfect crystal. A&R also allowed the diffusion of cations inside the contiguous phases, which is another aspect that has been observed experimentally. This diffusion is controlled by the randomization step during A&R and, in this particular case, it was minimized in order to avoid the formation of secondary phases or interdiffusion layers. Indeed, radial distribution functions of the cation-oxygen pairs demonstrated that no dopant clustering was formed during A&R. The small thickness of both layers, however, introduces strain effects derived from the lattice misfit between the materials. This strain provoked the compression of GDC and the expansion of YSZ, which could affect the oxygen conductivity. Oxygen ion diffusion as well as the polarization of the different ionic species will be investigated in future work.

| ASSOCIATED CONTENT
The electronic supporting information contains a short description of the theoretical background used in the METADISE program in order to obtain surfaces and grain boundaries, as well as a figure describing tilt and twisted grain boundaries.

ACKNOWLEDGMENTS
The authors acknowledge the Engineering and Physical Sciences Research Council (EPSRC) for financial support (grant references EP/K001329/1 and EP/K016288/1). 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 Program. NHdL thanks the Royal Society for an Industry Fellowship.

AUTHOR CONTRIBUTIONS
The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript.