John A.
Purton
*a,
Adam
Archer
b,
Neil L.
Allan
b and
David S. D.
Gunn
a
aScientific Computing Department, STFC, Daresbury Laboratory, Keckwick Lane, Warrington, WA4 4AD, UK. E-mail: john.purton@stfc.ac.uk
bSchool of Chemistry, Cantock's Close, University of Bristol, Bristol, BS8 1TS, UK
First published on 1st March 2016
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 Gd3+ segregates into domains. We have examined the local environment of the Gd3+ 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.
![]() | (1) |
For every two Gd3+ 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 often 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 significant 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 Drennan4 concluded the conducting properties are strongly influenced 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 fine 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 metal–oxygen 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 Gd3+ ions surround each oxygen vacancy and small such clusters are distributed in the ceria lattice at random, with no tendency of the Gd3+ ions to cluster together. Chen and Navrotsky9 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 Gd2O3 crystal structure (called C-type) in which the dopant metal atoms are 6-coordinate to oxygen and the fluorite (F-type) Raman mode gradually disappears when x increases beyond 0.2.10 The C-type Gd2O3 structure is closely related to that of fluorite requiring the removal of two oxygen atoms from every eight in an ordered fashion. A different ordering of vacancies gives rise to the pyrochlore structure – oxygen vacancies are all at next-nearest neighbour positions to Gd3+ and the coordination number of dopant cations is 8. The absence of clustering of Gd3+ is questioned11 by energy filtering transmission electron microscopy which has identified segregation of Gd3+ 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-filtering TEM (EFTEM) and selected-area electron diffraction (SAED)) have confirmed 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 finite 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 and trimers
and concluded that location of the Gd3+ 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 Gd3+ 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 next-nearest 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 nano-domains, 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 CeO2 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 CeO2 and Gd2O3; for x ≤ 0.25 the average structure is that of fluorite, 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 clarified. 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 configurations 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 limitations. Simulation cells containing many thousand atoms are feasible, and we can examine finite concentrations of dopants rather than solely clusters in the dilute limit. Temperature effects are included automatically. Thus here we present, to our knowledge, the first HMC study of grossly non-stoichiometric GDC.
The HMC approach we have developed is related to that used in the modelling of polymers and biomolecules.32–34 The technique has been applied previously to such problems as the enthalpies of mixing of binary oxides27,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 first of these is a short molecular dynamics (MD) simulation (140 steps and a timestep of 1.0 fs) in the microcanonical ensemble. The last configuration is accepted or rejected by comparing its energy with the energy of the starting configuration 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,
P(o→n) = exp{−β[Un − Uo]} |
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 configurations and takes explicit account both of ionic relaxation near impurity ions and thermal effects.
For consistency with previous calculations,16–22 the potential parameters used were those developed in ref. 13. A 20.0 Å cutoff for the potentials was employed.
A cation i (Gd or Ce) is selected from the structure. A set of even-order spherical harmonics40Qlm are calculated: for each vector r which connects this atom i to those within a specified primary cut-off distance,
Qlm(r) = Ylm(θ(r), ϕ(r)) | (2) |
The average value of Qlm over the Nb neighbours within the cut-off distance is thus
![]() | (3) |
We now average, for a given l, over all possible values of m to obtain Ql(i), which is rotationally invariant,
![]() | (4) |
To assign local structures surrounding individual atoms, we use the average of Ql(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 defines a secondary sphere and so we obtain finally Ql for each such sphere. Lechner and Dellago41 in their study of soft particle systems use a similar method but first find the mean of the Qlm. While this final averaging to obtain the Ql 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 Q2 (l = 2) and Q8 (l = 8).
![]() | (5) |
The ionic conductivity, σ, can then be obtained from the Nernst–Einstein relationship; note this is strictly only valid for dilute systems but in practice is often used also for solid solutions at the concentrations considered in this paper.
![]() | (6) |
We first present data for the pure end member oxides (CeO2 and Gd2O3) as this provides a platform to discuss the results for the CeO2/Gd2O3 mixtures. CeO2 has the fluorite structure (space group Fmm) in which all cation positions are symmetrically equivalent as are those of the anions. All the Ce–O and Ce–Ce (i.e. next nearest cation neighbours) distances are identical (Table 1).10 In contrast, in the C-type structure adopted by Gd2O3 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 fitted these to one or two Gaussians, as appropriate. In Fig. 1 we display the cation–cation RDFs for CeO2 and Gd2O3. As expected from the experimental data in Table 1 the M–M RDF for CeO2 has a single peak. This is in contrast to Gd2O3 where two peaks are observed (consistent with the bond lengths in Table 2) and this distinction can be used to interpret the structure of CeO2/Gd2O3 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.
Experiment | Calculated | |
---|---|---|
a | 5.405 | 5.442 |
Ce–O | 2.341 | 2.325 |
Ce–Ce | 3.822 | 3.824 |
O–O | 2.703 | 2.685 |
Experiment | Calculated | |
---|---|---|
a | 10.808 | 10.870 |
Gd(M1)–O | 2.349 | 2.291 |
Gd(M2)–O | 2 × 2.255 | 2.291 |
2 × 2.314 | ||
2 × 2.412 | ||
Gd–Gd | M1–M2 × 6 3.590 | 3.592/4.034 |
M1–M2 × 6 4.068 | ||
M2–M1 × 2 3.590 | ||
M2–M1 × 2 4.068 | ||
M2–M2 × 4 3.605 | ||
M2–M2 × 4 4.082 | ||
O–O | 2.976 |
![]() | ||
Fig. 1 Calculated cation–cation radial distribution functions (RDFs) for the nearest cation neighbours in (a) CeO2 and (b) Gd2O3. The red lines in figures (a) and (b) are fitted Gaussians. |
We now examine the change in lattice parameter as a function of composition. The solubility limit of Gd in CeO2 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 fluorite, while for 0.2 < x < 0.4 there is a two phase region of fluorite 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 fluorite 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 defined as the fluorite 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.
![]() | ||
Fig. 2 The calculated lattice parameter defined as the fluorite cell parameter, as explained in the text. The experimental data are from Chen and Navrotsky,9 Scavini et al.,23 Grover et al.,45 and Kennedy and Avdeev.44 |
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 co-workers.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 Gd3+ 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 Gd3+ 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. Moreover, experimental O–O bond lengths for GDC do the best of our knowledge have not been reported in the literature so no comparison is possible. However, there are O–O distances for Ce1−xYxO2−x/2 (YDC) from powder diffraction experiments47 demonstrating a similar trend and the data has been analysed using three distinct O–O bond lengths. The initial increase in the unit cell volume and subsequent decrease has been rationalised in terms of the relative importance of the O–O and cation–O bond lengths.23–25 The M–M distances, as discussed earlier, are a key marker of the local environment around the Ce4+ and Gd3+ ions. In Scavini et al.23 there are significant differences depending on the assumptions made to refine the data (i.e. whether the data was analysed using either the Rietveld or PDF methods). The Rietveld refinement for x = 0 can be reconciled with a single NNN distance consistent with the fluorite structure. As x increased the distance splits at xGd = 0.25 into short and long distances (Fig. 5). In contrast, in PDF analysis the longer M–M distance decreases slightly as xGd increases and is present at xGd < 0.25. Our results, and the EXAFS experiments in ref. 8 (also shown in Fig. 5), tend to support the PDF analysis. Indeed, for the short M–M distance our results are in excellent agreement with experiment and for the long M–M distance are intermediate between those of the PDF and EXAFS results. Even in the Ce–Ce RDF's (not shown) for xGd = 0.05 we observe a splitting of the peaks into two distinct distances. Our results do not support the formation of small pyrochlore-type clusters as suggested by Wang et al.20 since the Gd–Gd RDF for pyrochlore Gd2O3 does not show the splitting we see in Fig. 1b.
![]() | ||
Fig. 3 Ce–O distances as a function of composition for GdxCe1−xO2−x/2. The experimental data are from Ohashi et al.,5 Nakagawa et al.,6 and Yamazaki et al.7 |
![]() | ||
Fig. 4 Gd–O distances as a function of composition for GdxCe1−xO2−x/2. Experimental data are from Ohashi et al.,5 Nakagawa et al.,6 and Yamazaki et al.7 |
Composition x of GdxCe1−xO2−x/2 | O–O distance (Å) |
---|---|
0.00 | 2.685 |
0.05 | 2.695 |
0.10 | 2.701 |
0.15 | 2.721 |
0.20 | 2.735 |
0.30 | 2.753 |
0.40 | 2.750 |
0.50 | 2.796 |
0.60 | 2.818 |
0.70 | 2.850 |
0.80 | 2.898 |
0.85 | 2.921 |
0.90 | 2.940 |
0.95 | 2.957 |
1.00 | 2.976 |
![]() | ||
Fig. 5 The M–M distances as a function of composition. The PDF and Rietveld analysis are from ref. 23 and the EXAFS data from ref. 8. |
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 Gd3+ defects.48 Thus, in order to demonstrate the clustering of Gd3+ ions we have removed all the Ce4+ and O2− ions in Fig. 6 to display only the Gd3+ 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 Gd3+ ions. In addition, we have calculated both the average ratio of isolated to the total number of Gd3+ ions (Fig. 7a) and total number (Fig. 7b) of isolated Gd3+ ions (i.e. an isolated Gd3+ ion is where there are no Gd3+ ions as nearest cation neighbours). At x = 0.05 approximately half of the Gd3+ ions are isolated, the remaining Gd3+ ions have at least one Gd3+ 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–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 combinations13 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 Gd3+ ions increases (Fig. 6b and 7), the fraction of isolated Gd3+ ions decreases. Further increases in Gd content cause a dramatic decrease in both the average ratio and absolute total of isolated Gd3+ ions. Indeed, at x = 0.2 there appears to be the formation of a network consisting largely of Gd3+ ions and by x = 0.3 few isolated Gd3+ ions remain.
To obtain more information on local structure in the solid solution, Q2 and Q8 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 CeO2, bulk cubic Gd2O3 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 for the compositions x = 0.1, x = 0.2, x = 0.4 and x = 0.8; other compositions are consistent with these results. The signatures for each composition are distinct and there is an evident shift from the CeO2 end-member towards the other end-member, Gd2O3, as the Gd content increases. The change in peak position with composition is not linear. The width of the signatures reflects the variation in local environments for both Gd and Ce, which are reflected in the O–Gd–O and O–Ce–O angles respectively. It is evident from Fig. 8 that CeO2 and Gd2O3 produce distinct features in the Steinhardt order parameters that reflect the average local environment of the Ce4+ and Gd3+ cations. The order parameters are dissimilar to those of other oxide structures such as the pyrochlore Gd2Ce2O7 (for which the Q2 and Q8 peaks are narrower and at 0.15 and 0.40 respectively) and corresponding defective fluorite and the signatures can be employed as an indication of the local environment of the cations in the GdxCe1−xO2−x/2 mixtures. Starting from the CeO2 end member, as Gd/O dopants/vacancies are created the order parameters of both Ce4+ and Gd3+ are positioned close to that of the CeO2 end-member, reflecting local environments similar to that in CeO2. The fluorite framework dominates. The further introduction of Gd causes a local distortion of both cations. The widths of the signatures for Gd3+ are broader than for Ce4+ indicating that the range of angles around the Gd3+ ions is much greater than that around the Ce4+ ions, reflecting the larger Ce–O bond strength. As more Gd3+ ions are added the probability distributions become progressively more similar to that of Gd2O3. In addition, the width of the Gd3+ peak decreases, whilst that of the Ce4+ gradually increases. The structure becomes more similar to that of Gd2O3.
![]() | ||
Fig. 8 Probability distributions of the order parameters obtained from the hybrid Monte Carlo simulations. (a) Q2 for Ce, (b) Q2 for Gd, (c) Q8 for Ce and (d) Q8 for Gd. |
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: (A) cells in which the Gd, and independently the oxygen vacancies, are distributed at random; (B) cells in which all the Gd ions and oxygen vacancies are distributed around the centre of the cell, forming one discrete nano-domain, or nano-crystallite of Gd2O3. Two different such examples for (B) were examined, in which the interfaces of the nanocrystallite were {100}, and {111} respectively. Fig. 9 displays the Q2 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 artificially constructed cells. For the (B) cells the signatures are further away from the Gd2O3 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. Q8 plots show the same variations. Overall these results suggest the formation of smaller nanodomains of Gd2O3 in the HMC simulations with a greater uniformity of Gd environments.
Visualisation of the positions of cations obtained from the simulations reveal that the Gd3+ ions start to form domains even at low concentrations (xGd ≈ 0.1). Moreover, we have demonstrated that Gd3+ rich domains are established when xGd reaches approximately 0.2. This is in contrast to the solubility limits of Gd3+ in CeO2 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 Gd3+ rich network reduces ionic conductivity. Moreover, we propose that the growth in the network of Gd3+ 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 Gd2O3 nanodomains restrict the diffusion paths and decrease the mobility of oxygen vacancies. Fig. 6 shows that by x = 0.3 few isolated Gd3+ 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 Gd2O3-like 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 influenced 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.
Footnote |
† 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. |
This journal is © The Royal Society of Chemistry 2016 |