Growth of nano-domains in Gd–CeO2 mixtures: hybrid Monte Carlo simulations

Hybrid Monte Carlo (HMC) simulations are used to study the growth of Gd-rich domains in Gd doped CeO2, and we probe the conductivity of the resulting and other configurations by molecular dynamics. Previous work has been restricted to the dilute defect limit, assumptions of particular cluster formation, and neglect of all temperature effects. Our methods suffer none of these restrictions. Even at low concentrations Gd segregates into domains. We have examined the local environment of the Gd ions using radial distribution functions and Steinhardt order parameters. The observed structure is consistent with the formation of cubic C-type (Gd2O3) domains, rather than the monoclinic B-type or pyrochlore clusters which have been suggested previously. In addition, previous detailed pair distribution function analysis of the solid solution has indicated different local cation environments from those from a Rietveld analysis – overall our results support the former analysis rather than the latter. At the elevated temperatures (1000 K) of the simulations there is no particular preference for vacancy and dopant cations to be located at second neighbour sites, an issue long discussed for this and similar systems. Both calculated and experimental conductivities show a similar variation with composition, passing through a maximum with increasing Gd concentration. The conductivities of the configurations generated in the hybrid Monte Carlo simulations are lower than those of configurations generated independently in which the Gd ions are distributed at random. The HMC thermally generated Gd nano-domains capture oxygen vacancies, reduce oxygen vacancy mobility and block diffusion paths.


Introduction
CeO 2 is a technologically important material, abundant and with rich defect-and nano-chemistries. 1 In particular gadolinium doped ceria (GDC), Gd x Ce 1Àx O 2Àx/2 , has been proposed as an intermediate temperature solid oxide fuel cell (SOFC) electrolyte, 2 and consequently has been the subject of numerous experimental and theoretical studies. The efficiency of a SOFC is linked to the ionic conductivity and the conductivity of CeO 2 can be enhanced by doping with a cation with a different charge to that of Ce. Here Gd 3+ replaces Ce 4+ : For every two Gd 3+ dopants introduced, an oxygen vacancy is created and eqn (1) implies that the ionic conductivity should increase steadily with dopant concentration and thus proportional to the concentration of oxygen vacancies. However the actual behaviour is considerably more complex and the experimental conductivity of GDC is not a simple function of the dopant concentration, but passes through a maximum at approximately 10-20% Gd, above which the conductivity starts to decrease. 2 The conductivity is also a function of the thermal history of the sample. Zhang et al., for instance, have compared the conductivity of sintered pellets and samples aged at 1273 K for 3 to 8 days. 3 For compositions x $ 0.2 they observed a marked decrease in conductivity for the aged samples and attributed this to the formation of micro-domains. The experimental conductivity is oen further complicated by the presence of grain boundaries; the dc conductivity will also include contributions from grain boundaries, and production of an optimum electrolyte has to involve optimisation of the microstructure.
This paper concentrates on the properties of the bulk material. The variation of the conductivity with Gd concentration has long been thought to be due to interactions between the dopant cations and oxygen vacancies and the formation of associated clusters. This is a challenging problem for both experiment and theory and has received signicant and continuing attention due to the importance of understanding the mechanisms of conductivity and potential improvement in the performance of SOFC. High resolution transmission electron microscopy (HRTEM) carried out by Mori and Drennan 4 concluded the conducting properties are strongly inuenced by the micro-domain size in the grain; the larger the microdomains the lower the conductivity and the authors suggested the domains may have a distorted pyrochlore structure. From extended X-ray adsorption ne structure (EXAFS) spectra, probing the local structure around cations, Ohashi et al. 5 suggested the association of two Gd ions and one oxygen vacancy in heavily doped solutions with vacancies located around both Gd and Ce. XAFS studies in ref. 6 and 7 concluded that metaloxygen bond lengths increase with decreasing Gd content. Deguchi et al. 8 used EXAFS to examine the coordination of dopant cations and suggested that 1-2 Gd 3+ ions surround each oxygen vacancy and small such clusters are distributed in the ceria lattice at random, with no tendency of the Gd 3+ ions to cluster together. Chen and Navrotsky 9 combined calorimetry and X-ray diffraction, and in line with Deguchi et al. 8 propose the same neutral trimers as the main associated defects and speculate on the importance of local site distortion. In contrast Banerji et al. interpreted their experimental Raman spectra in terms of the formation of regions with the Gd 2 O 3 crystal structure (called C-type) in which the dopant metal atoms are 6coordinate to oxygen and the uorite (F-type) Raman mode gradually disappears when x increases beyond 0.2. 10 The C-type Gd 2 O 3 structure is closely related to that of uorite requiring the removal of two oxygen atoms from every eight in an ordered fashion. A different ordering of vacancies gives rise to the pyrochlore structureoxygen vacancies are all at next-nearest neighbour positions to Gd 3+ and the coordination number of dopant cations is 8. The absence of clustering of Gd 3+ is questioned 11 by energy ltering transmission electron microscopy which has identied segregation of Gd 3+ to nano-sized domains even for x ¼ 0.1. The size of these domains increased with increasing x and the analysis of Ye et al. 11 who also revise the earlier conclusions of ref. 4, suggests a C-type rather than pyrochlore structure for the domains. The measured electrical conductivity led the authors to speculate that the conductivity is inversely proportional to their size. 11 More recent microscopy (HRTEM, energy-ltering TEM (EFTEM) and selected-area electron diffraction (SAED)) have conrmed the existence of defect clusters with ordered structures in GDC. 12 In parallel there has been substantial computational work to try and resolve some of the issues raised by these different experimental conclusions. For methodological reasons to which we return below direct contact between experiment and theory is not straightforward. Thus, early studies have concentrated on point defect calculations and energy minimisation using interatomic potentials for particular clusters in the dilute limit rather than solid solutions with nite dopant concentrations. These studies have also been restricted to the static limit, i.e., 0 K and neglecting vibrational contributions to the energy. Minervini et al. 13 using such methods examined different arrangements of isolated dimers ðGd 0 Ce : V cc O Þ c and trimers ð2Gd 0 Ce : V cc O Þ Â and concluded that location of the Gd 3+ at nearest neighbour † (NN) or next-nearest neighbour (NNN) positions to the oxygen vacancies to were very similar in energy. More recent ab initio DFT+U calculations however have indicated the nearest neighbour arrangement is lower in energy. 14, 15 Ye et al. have used the same computational approach as Minervi et al., 13 with a different set of interatomic potentials, for larger, but still isolated, defect clusters with up to six vacancies and associated dopant cations. 16 In their most stable clusters containing more than three oxygen vacancies, each vacancy has two Gd 3+ cations in the next-nearest neighbour site common to both, in a C-type consistent structure. However, similar studies by Li et al. 17,18 favour nearest rather than nextnearest neighbours. Lower energy clusters than those in ref. 16 have been reported subsequently by Li et al. 19 Again using the same computational methods, Wang et al. 20 suggest that dopant concentration limits the size of the nanodomains, such that at lower Gd concentrations GDC favours small pyrochlore-type clusters but C-type nano-domains at higher concentration and suggest the maximum in the conductivity vs. concentration plot is linked to the formation of the latter. In a related study Wang et al. 21 conclude from calculated (dilute limit) defect energies that monoclinic B-type extended defect structures within the CeO 2 lattice are less stable than C-type. Wang and Cormack have more recently also investigated the relationship between bulk strain and defect structure in GDC and proposed that the next nearest neighbour type associations are preferred under zero or tensile applied strains, while nearest neighbour type defect associations become more favourable under compression. 22 Evidence for the presence of extended defects locally has come from an experimental Pair Distribution Function (PDF) Analysis by Scavini et al. 23 However, as Fig. 4 of ref. 18 shows, there are major differences in the variation of the interatomic metal-metal and metal-oxygen interatomic distances determined by Rietveld analysis, PDF and EXAFS. We return to this discrepancy later in this paper. Most recently Scavini et al. 24 conclude that on the subnanometre scale all GDC solid solutions can be described as biphasic containing 'droplets' of CeO 2 and Gd 2 O 3 ; for x # 0.25 the average structure is that of uorite, which turns into C-type at larger x and they view these changes as a percolation-driven phase transition.
Thus in spite of a large number of experimental and computational studies, the nature and role of extended defect clusters in GDC and how these vary with composition are far from being claried. A study providing a complete description of both local and micro-structure is required. Moreover, the atomistic simulations described above employed energy minimisation in the static limit (i.e., at 0 K, ignoring vibrational terms). Absolute and relative defect binding energies can strongly vary with temperature; Raman spectroscopy of GDC at elevated temperatures indicates that oxygen vacancies dissociate from clusters stable at lower temperatures. 25 Indeed, the number of congurations and the size of the defect clusters explored in previous computational studies are limited and this in turn has restricted the contact between theory and experiment.
In this paper we report hybrid Monte Carlo (HMC) simulations, which combine Monte Carlo and molecular dynamics techniques in order to explore phase space efficiently. These are an appropriate method to remove almost all of these † We use, following the literature, NN throughout this paper to refer to a cation at a nearest neighbour position to an oxygen vacancy and NNN to refer to a cation at a next-nearest (cation) position to an oxygen vacancy.
limitations. Simulation cells containing many thousand atoms are feasible, and we can examine nite concentrations of dopants rather than solely clusters in the dilute limit. Temperature effects are included automatically. Thus here we present, to our knowledge, the rst HMC study of grossly nonstoichiometric GDC.

Monte Carlo simulations
The properties of solid state materials, especially, ionic compounds, have been investigated using either supercell or point defect calculations. 26 However, these methods are not readily extended to disordered systems containing a nite impurity or defect content, thereby restricting simulations to end-member compounds and excluding studies of many naturally occurring minerals and ceramics of industrial importance. Standard Monte Carlo (MC) and molecular dynamics (MD) is unable to overcome kinetic energy barriers to sufficiently sample all necessary congurations. In previous papers we have described MC methods that exchange cation positions in order to sample multiple congurations and calculate the thermodynamic and solubility limits of ionic materials. [27][28][29][30][31] However, even with MC swaps the sampling of congurations will be very poor due to the different size and charge of the cations. It is necessary to resort to a different strategy, i.e. some form of relaxation must be provided. We also demonstrated, in a previous paper, that a simple relaxation (energy minimisation) is not sufficient and vibrational effects must be included within the calculation. 27 Unfortunately these effects were not included in previous calculations.
The HMC approach we have developed is related to that used in the modelling of polymers and biomolecules. [32][33][34] The technique has been applied previously to such problems as the enthalpies of mixing of binary oxides 27,29 and non-convergent ordering of cations in olivine and phase transitions. 19 During one HMC cycle, one of three options is chosen at random, with equal probability. The rst of these is a short molecular dynamics (MD) simulation (140 steps and a timestep of 1.0 fs) in the microcanonical ensemble. The last conguration is accepted or rejected by comparing its energy with the energy of the starting conguration and using the standard Metropolis algorithm. 35 The Metropolis approach requires that a trial move from the original state (o) to a new state (n) is accepted with the probability, where U n and U o are the energies of the new and original states and b is 1/kT. At the start of each MD run, velocities are chosen anew at random from a Maxwellian distribution. In the second, which is only applicable to the solid solution, a short MD run follows a random exchange of Ce and Gd ions. Again, the difference in energy between the previous conguration and that immediately aer the MD simulation is used in the Metropolis algorithm. Since the MD technique is used to update the positions of the ions it not possible to track the location of the vacancies and it was not feasible to swap the position of a vacancy with an oxygen ion. However, in GDC the oxygen ions are mobile at 1000 K and can move to the most energetically favourable position. The third option is a random change of the volume/shape of the box, which again is accepted or rejected using the Metropolis algorithm. Enthalpy and structural data were averaged over a period of 150 000 cycles, prior to which an equilibration period of 150 000 cycles was undertaken.
We stress that our method does not involve the use of an approximate parameterised Hamiltonian. Not only does parameterisation of, for example, an Ising-type Hamiltonian become increasingly difficult beyond binary or pseudobinary mixtures, but it can average out local effects due to ion clustering and association, and such methods cannot readily be extended to include the effects of lattice vibrations and pressure. The HMC technique permits an efficient sampling of different congurations and takes explicit account both of ionic relaxation near impurity ions and thermal effects.

Order parameters
To examine the nature of the local order more closely we have evaluated averaged local (Steinhardt) bond order parameters. 36 These are oen used to distinguish different crystalline structures and polytopes, (e.g. ref. 37), to characterise non-stoichiometry at interfaces 38 and very recently the disorder induced by radiation damage and subsequent healing in ceramics. 39 These parameters are translationally and rotationally invariant, i.e., independent of the reference frame specifying the crystal structure.
A cation i (Gd or Ce) is selected from the structure. A set of even-order spherical harmonics 40 Q lm are calculated: for each vector r which connects this atom i to those within a specied primary cut-off distance, q and f are polar angles which dene the orientation of r. Here we set the primary cut-off to 3.2Å, which restricts contributions to Q lm from only nearest oxygen neighbours. The average value of Q lm over the N b neighbours within the cut-off distance is thus We now average, for a given l, over all possible values of m to obtain Q l (i), which is rotationally invariant, To assign local structures surrounding individual atoms, we use the average of Q l (i) over all atoms i located within 8Å, the secondary cut-off, of one of the atoms in the simulation cell. The region around each atom i denes a secondary sphere and so we obtain nally Q l for each such sphere. Lechner and Dellago 41 in their study of so particle systems use a similar method but rst nd the mean of the Q lm . While this nal averaging to obtain the Q l invariably reduces the spatial resolution, it makes it easier in practice to distinguish different structures, and tests showed that 8Å was the optimum secondary radius to achieve this aim. In this paper we use Q 2 (l ¼ 2) and Q 8 (l ¼ 8).

Molecular dynamics
In order to compare with experimental conductivity measurements calculation of the ionic diffusion is required. This cannot be achieved using the HMC technique and we have chosen to use molecular dynamics simulations. Atomic congurations obtained from the HMC were used as starting points for the simulations. The MD simulations, with the DL_POLY package, 42 were undertaken in the isothermal-isobaric ensemble using Nosé-Hoover dynamics to control the temperature (1000 K) and pressure (1 atmosphere). The oxygen diffusion constant, D, was calculated from the mean squared displacements following the procedure in ref. 43: The ionic conductivity, s, can then be obtained from the Nernst-Einstein relationship; note this is strictly only valid for dilute systems but in practice is oen used also for solid solutions at the concentrations considered in this paper.
where z i e is the charge of species i, and c is the concentration of defectshere oxygen vacancies, which we assume is the only species responsible for the conductivity.

Hybrid Monte Carlo
The   (Table 1). 10 In contrast, in the C-type structure adopted by Gd 2 O 3 there are two distinct cation and two distinct anion sites (M1 and M2 for the cations, O1 and O2 for the anions, respectively). The M1 site has six Gd-O bonds of the same length, but the M2 is environment is less symmetric, giving rise to a range of bond lengths (Table 2) 10,44,45 and this distortion also gives rise to a range of Gd-Gd distances. In our simulations the automatic calculation of bond lengths is straightforward and an appropriate vehicle for their presentation is the radial distribution function (RDF). We have calculated the RDFs for HMC simulations and to obtain average bond lengths we have tted these to one or two Gaussians, as appropriate. In Fig. 1 we display the cation-cation RDFs for CeO 2 and Gd 2 O 3 . As expected from the experimental data in Table 1 the M-M RDF for CeO 2 has a single peak. This is in contrast to Gd 2 O 3 where two peaks are observed (consistent with the bond lengths in Table 2) and this distinction can be used to interpret the structure of CeO 2 /Gd 2 O 3 solid solutions as discussed below. The lattice parameter and bond lengths determined from X-ray diffraction are compared with data determined using HMC simulations in Tables 1 and 2. Our simulations have a slightly greater lattice parameter than experiment (0.7% difference), which, in part, can be attributed to the higher temperature of the simulation than experiment.
We now examine the change in lattice parameter as a function of composition. The solubility limit of Gd in CeO 2 has been investigated using X-ray diffraction and reported by Bevan et al. to exceed 40%. 46 This has been disputed by Chen and Navrotsky who observed additional peaks in their diffraction patterns due to C-type phase at 39% and the lattice parameter varies linearly with Gd content only up to 30.5%. 9 Chen and Navrotsky interpret these diffraction data to infer that at x # 0.2 the structure is uorite, while for 0.2 < x < 0.4 there is a two phase region of uorite and C-type and for x > 0.4 a C-type solid solution. Comparison of our calculations and experiment is difficult; not only were our calculations run at 1000 K but the data can be interpreted assuming several possible models. Thus in ref. 9 the authors have presented data that has been analysed within the above constraints and so for x > 0.4 split the uorite and C-type lattice parameters. Indeed, Scavini et al. concluded that diffraction data could not be correctly interpreted using average structural models but required a biphasic model even for x ¼ 0.125. 23 Unfortunately, making a direct comparison with the biphasic model interpretation of the experimental data is not possible as we are only able to determine the lattice parameters from the simulation cell vectors and cannot identify individual components in the same manner as XRD which operates over a much longer length scale. For the sake of comparison we compare our calculated values with experimental data dened as the uorite cell parameter (i.e. half the C-type parameter) (Fig. 2). The calculated values are greater than the experimental values by less than 1%. Both experimental and calculated values exhibit a strong positive deviation from Vegard's law, i.e. from ideal solid solution behaviour. The calculated values form a slightly asymmetric curve with a maximum at x ¼ 0.6. We turn to analyse trends in bond lengths and compare those determined using EXAFS and in a detailed pair distribution function (PDF) and Rietveld analysis by Scavini and coworkers. 23 We commence with the cation-oxygen (i.e. nearest neighbour) bond lengths. The variation of the Ce-O (Fig. 3) and of the Gd-O (Fig. 4) bond lengths reproduce the observed EXAFS trends although underestimated by the current potential model. For the end-member compounds the difference in the average bond lengths for Ce-O and Gd-O are 0.5% and 2% respectively. The contraction in cation-oxygen bond lengths with Gd content is in contrast to the lattice parameter in which increases with Gd 3+ concentration. The average O-O bond distances obtained from our simulations is presented in Table 3 and exhibit a gradual increase as the number of Gd 3+ ions increases. We stress that the relatively high temperature of the simulations and the distortion of octahedra associated with the formation of interfaces produces a wide range of O-O distances that can only be interpreted in terms of a simple average.     23 there are signicant differences depending on the assumptions made to rene the data (i.e. whether the data was analysed using either the Rietveld or PDF methods). The Rietveld renement for x ¼ 0 can be reconciled with a single NNN distance consistent with the uorite structure. As x increased the distance splits at x Gd ¼ 0.25 into short and long distances (Fig. 5). In contrast, in PDF analysis the longer M-M distance decreases slightly as x Gd increases and is present at Although we cannot make direct comparison with the XRD lattice parameters, we are able to examine snapshots of the structure and in particular the distribution of the cations. Several authors have postulated that the ionic conductivity is related to the trapping of vacancies around clusters of Gd 3+ defects. 48 Thus, in order to demonstrate the clustering of Gd 3+ ions we have removed all the Ce 4+ and O 2À ions in Fig. 6 to display only the Gd 3+ ions and where any two Gd are nearest (cation) neighbours to each other a "bond" has been inserted (this "bond" is only for graphical purposes and there no other inference as to physical interaction)we call these associated Gd 3+ ions. In addition, we have calculated both the average ratio of isolated to the total number of Gd 3+ ions (Fig. 7a) and total number (Fig. 7b) of isolated Gd 3+ ions (i.e. an isolated Gd 3+ ion   is where there are no Gd 3+ ions as nearest cation neighbours). At x ¼ 0.05 approximately half of the Gd 3+ ions are isolated, the remaining Gd 3+ ions have at least one Gd 3+ ion as a nearest (cation) neighbour (Fig. 6a). We have already referred to a number of computational studies which have demonstrated the favourable formation of vacancy-dopant clusters in doped ceria. [16][17][18][19][20][21][22] These calculations were performed at 0 K and for the potential model employed in this study the difference in defect binding energy between these combinations 13 is relatively low for Gd, as discussed above. In contrast our calculations were run at 1000 K and so explore a wider range of defect combinations without any a priori assumptions. As x increases to 0.10 while the total number of isolated Gd 3+ ions increases (Fig. 6b and 7), the fraction of isolated Gd 3+ ions decreases. Further increases in Gd content cause a dramatic decrease in both the average ratio and absolute total of isolated Gd 3+ ions. Indeed, at x ¼ 0.2 there appears to be the formation of a network consisting largely of Gd 3+ ions and by x ¼ 0.3 few isolated Gd 3+ ions remain.
To obtain more information on local structure in the solid solution, Q 2 and Q 8 Steinhardt signatures were calculated with Gd at the centre of the primary sphere and separately also with Ce from 50-100 snapshots from the equilibrated hybrid Monte Carlo runs at 1000 K on bulk CeO 2 , bulk cubic Gd 2 O 3 and on each composition. All individual curves are normalised such that the area under each curve is one. For Gd and Ce these are collected together in Fig. 8  To investigate further the sensitivity of the order parameter analysis to different arrangements of the Gd and the oxygen vacancies we have also analysed separately the results of optimised supercells of the same composition and size but which were not generated in the HMC runs. These included:  Fig. 9 displays the Q 2 and Steinhardt signature for all these cases. None of these match the signatures from the HMC; it is rather hard to judge but (A) is closer than (B) to the HMC. Almost all the signatures are wider than those in Fig. 8 from the HMC, showing a greater spread of angles and more distorted environments in the articially constructed cells. For the (B) cells the signatures are further away from the Gd 2 O 3 signature than those of the HMC plots, even with increasing Gd concentration, which must be a consequence of the different boundaries (interfaces) separating Gd and Ce rich regions. Q 8

Conductivity
Previous experimental and theoretical research, as discussed in the introduction, has identied a complex relationship between the atomic conguration and the ionic conductivity and the formation of Gd-rich domains reduces the ionic conductivity. In ref. 3 the conductivity was measured in freshly sintered (1873 K for 5 hours) material and in samples that had been aged for up to 8 days (heating at 1273 K). The ratio of the conductivity of the aged and sintered samples demonstrated that for x < 0.2 the aged samples had the greater conductivity whilst for x > 0.2 the aged samples had signicant lowering of the conductivity. This observation was explained by the precipitation and coarsening of these domains during the ageing process. 3 To test this hypothesis and the inuence of Gd-rich domains on the ionic conductivity, we have undertaken MD simulations, as described in Section 2.3, using different cation congurations as starting points for the simulations. We have assumed that the sintered sample can be represented by a random conguration of cations (and vacancies) and the aged samples by congurations obtained from the HMC simulations (representative congurations were taken from the later stages of the calculations). The ratio of the oxygen ion conductivity obtained from the random, s random , and HMC, s HMC , congurations is presented in Fig. 10. Our results are in good agreement with the experimental results of ref. 3 and the calculated conductivity of the HMC congurations is greatly reduced for x $ 0.15. This is the composition   10 The ratio of ionic conductivity obtained from HMC starting configurations, s HMC , and those obtained from random configurations of cations, s random . The ionic conductivity is assumed to be entirely from oxygen ion diffusion.
at which we observe the formation of a network of Gd-rich ions. Despite the signicant approximations in our calculations, our results suggest that the growth of Gd-rich domains have signicant impact on the conductivity of GDC.

Conclusions
Hybrid Monte Carlo calculations have been performed to examine the thermodynamic equilibrium properties of Gd doped ceria. In particular, HMC simulations are capable of providing efficient sampling of congurations where large differences in size of the cations/anions. Unlike previous calculations 13, [17][18][19][20][21][22] our work includes the effect of temperature on defect associations and the HMC technique allows exploration of a large number of congurations with different cation arrangements. We make no assumptions as to the formation of any particular clusters, and temperature and effects are included automatically and we are not restricted to the dilute limit. It is worth stressing our HMC simulations do not take into account any kinetic factors, which are undoubtedly important in the fabrication of experimental samples and their behaviour over long time scales.
Visualisation of the positions of cations obtained from the simulations reveal that the Gd 3+ ions start to form domains even at low concentrations (x Gd z 0.1). Moreover, we have demonstrated that Gd 3+ rich domains are established when x Gd reaches approximately 0.2. This is in contrast to the solubility limits of Gd 3+ in CeO 2 reported by Bevan et al. where solid solution systems greater than 40% Gd were observed. 46 We stress that the hybrid Monte Carlo technique explores equilibrium thermodynamic properties. This is in contrast to many experiments are unlikely to have been carried out under equilibrium because of the very long annealing times necessary.
MD simulations have been used to calculate the ionic conductivity on the representative structures obtained from the HMC. The results of these calculations and those in which the positions of the ions are randomised support the hypothesis that a Gd 3+ rich network reduces ionic conductivity. Moreover, we propose that the growth in the network of Gd 3+ ions is the main driver in the reduction of the conductivity of GDC. Overall our results are consistent with percolation theory arguments; 49 as the Gd concentration increases the Gd 2 O 3 nanodomains restrict the diffusion paths and decrease the mobility of oxygen vacancies. Fig. 6 shows that by x ¼ 0.3 few isolated Gd 3+ ions remain and this is close to the site percolation threshold of 0.311 for a cubic lattice, and where Scavini et al. suggest a percolation-driven phase transition takes place. 24 We have analysed both the RDF's and Steinhardt order parameters to examine the local environment around the cations, and both are consistent with the formation of Gd 2 O 3like regions. The calculated cation-cation RDF's lend support to the interpretation of powder diffraction data using PDF's rather than the more standard Rietveld analysis of the data. 10 The calculated RDF's are also in reasonable agreement with data obtained from EXAFS experiments. 7,8 The shorter than expected cation-oxygen bond length suggests that these potentials are too strong and may have inuenced previous simulations concerned with defect binding. Further work will improve the efficiency of the technique so that we can study much larger simulation cells to observe growth of domain structure in materials.