# Direct interpretation of the X-ray and neutron three-dimensional difference pair distribution functions (3D-ΔPDFs) of yttria-stabilized zirconia

Ella Mara Schmidt, Reinhard B. Neder, James D. Martin, Arianna Minelli, Marie-Hélène Lemée, Andrew L. Goodwin

► To cite this version:

Ella Mara Schmidt, Reinhard B. Neder, James D. Martin, Arianna Minelli, Marie-Hélène Lemée, et al.. Direct interpretation of the X-ray and neutron three-dimensional difference pair distribution functions (3D-ΔPDFs) of yttria-stabilized zirconia. Acta Crystallographica Section B: Structural Science, Crystal Engineering and Materials [2014-..], 2023, 79 (2), pp.138-147. 10.1107/S205252062300121X. hal-05442034

HAL Id: hal-05442034

https://hal.science/hal-05442034v1

Submitted on 13 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

Acta Cryst.

B

STRUCTURAL SCIENCE
CRYSTAL ENGINEERING
MATERIALS

ISSN 2052-5206

Received 21 December 2022

Accepted 9 February 2023

Edited by R. Černý, University of Geneva,
Switzerland

**Keywords:** diffuse scattering; cubic zirconia; three-dimensional difference pair distribution function (3D-ΔPDF).

**Supporting information:** this article has
supporting information at journals.iucr.org/b

# Direct interpretation of the X-ray and neutron three-dimensional difference pair distribution functions (3D-ΔPDFs) of yttria-stabilized zirconia

Ella Mara Schmidt,<sup>a,b*</sup> Reinhard B. Neder,<sup>c</sup> James D. Martin,<sup>d</sup> Arianna Minelli,<sup>b</sup> Marie-Hélène Lemée<sup>e</sup> and Andrew L. Goodwin<sup>b</sup>

Three-dimensional difference pair distribution functions (3D-ΔPDFs) from X-ray and neutron diffraction experiments are reported for yttria-stabilized zirconia (Zr<sub>0.82</sub>Y<sub>0.18</sub>O<sub>1.91</sub>). A quantitative analysis of the signatures in the three-dimensional difference pair distribution functions is used to establish that oxygen ions neighbouring a vacancy shift by 0.525 (5) Å along ⟨1, 0, 0⟩ towards the vacancy while metal ions neighbouring a vacancy shift by 0.465 (2) Å along ⟨1, 1, 1⟩ away from the vacancy. The neutron 3D-ΔPDF shows a tendency for vacancies to cluster along (1/2, 1/2, 1/2), which results in sixfold coordinated metal ions.

## 1. Introduction

The technological importance of zirconia (ZrO₂) is undoubtedly and this is underlined by its numerous applications in the ceramic industry. At ambient pressure and temperature ZrO₂ is monoclinic (P2₁/c) but upon increasing temperature it transforms to a tetragonal phase (P4₂/nmc, at approximately 1440 K) and a cubic phase (Fm$\bar{3}$m, at approximately 2640 K) (Boysen et al., 1991; Bondars et al., 1995). The addition of aliovalent oxides (CaO, MgO or Y₂O₃) allows the cubic phase to be stabilized at ambient conditions and results in a material that is widely applied because of its strength, high refractive index and thermal-shock resistance (Stapper et al., 1999). To maintain the overall charge balance in cubic stabilized zirconia (CSZ), oxygen vacancies are introduced into the structure which makes CSZ an ion conductor used in solid oxygen fuel cells (Tsampas et al., 2015), oxygen sensors (Schindler et al., 1989) and oxygen pumps (Pham & Glass, 1998).

The high-temperature cubic polymorph of pure ZrO₂ adopts the fluorite structure, where each Zr⁴⁺ is in regular eightfold coordination. Upon cooling, the local coordination environment in pure ZrO₂ is distorted and Zr⁴⁺ is sevenfold coordinated. On average, CSZ still adapts the fluorite structure but the introduction of oxygen vacancies reduces the average coordination number to 8 - x, where x depends on the dopant ion charge and concentration. The resulting complex defect structure of CSZ is a long-standing problem that has been addressed by numerous computational and experimental means. The nature of the most localized interactions is widely agreed upon in the literature and reproduced in Fig. 1 (Frey et al., 2005; Khan et al., 1998; Fèvre et al., 2005): oxygen ions neighbouring a vacancy relax along ⟨1, 0, 0⟩ towards this vacancy. Electrostatic effects and atomic sizes

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

OPEN

<figure><img src="image_5.png" /></figure>

ACCESS

Published under a CC BY 4.0 licence

<footnote><sup>a</sup>Faculty of Geosciences, MARUM and MAPEX, University of Bremen, Bremen, Germany, <sup>b</sup>Inorganic Chemistry</footnote>

<footnote>Laboratory, University of Oxford, Oxford, United Kingdom, cInstitute of Condensed Matter Physics, Friedrich-Alexander</footnote>

<footnote>University, Erlangen, Germany, dDepartment of Chemistry, North Carolina State University, Rayleigh, USA, and</footnote>

<footnote>eInstitut Laue-Langevin, Grenoble, France. *Correspondence e-mail: ella.schmidt@uni-bremen.de</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 1**
Schematics of local relaxations in YSZ as reported in the literature (Frey *et al.*, 2005; Khan *et al.*, 1998; Fèvre *et al.*, 2005). Oxygen ions (red) relax along (1, 0, 0) towards a neighbouring vacancy, next neighbour (NN) Zr ions relax along (1, 1, 1) away from the vacancy, while next nearest neighbour (NNN) dopant ions do not show significant relaxations.</figcaption></figure>

need to be taken into account when evaluating which type of metal ions neighbours a vacancy. Whether the oxygen vacancies are located preferably as a nearest neighbour (NN, i.e. in the coordination of the dopant metal ion) or a next nearest neighbour (NNN) depends on the concentration, size and charge of the dopant ion (Khan *et al.*, 1998). For yttria-stabilized zirconia (YSZ), simulations show that NNN vacancies are preferred (Bogicevic & Wolverton, 2003; Khan *et al.*, 1998). The $Zr^{4+}$ ions next to the vacancy relax along (1, 1, 1) away from the vacancy (Frey *et al.*, 2005; Khan *et al.*, 1998; Fèvre *et al.*, 2005). The exact magnitude of the relaxations and correlated further neighbour displacements depend on type and concentration of the dopant ion and only a limited agreement is reported in literature [see Frey *et al.* (2005) for a review].

In our contribution we focus on the defect structure in YSZ because of its technological importance as an oxygen ion conductor. We use three-dimensional difference pair distribution functions (3D-ΔPDFs) from neutron and X-ray diffraction experiments to quantify local-order principles. The material under investigation is by no means a newly studied material because YSZ was first described in 1951 by Hund (1951). Previous theoretical and experimental investigations include first-principle calculations (Stapper *et al.*, 1999), molecular dynamics simulations (Bogicevic & Wolverton, 2003; Fabris *et al.*, 2002), Bragg data analysis (Kaiser-Bischoff *et al.*, 2005; Ishizawa *et al.*, 1999; Morinaga *et al.*, 1979), extended X-ray absorption fine structure (EXAFS) (Ishizawa *et al.*, 1999; Catlow *et al.*, 1986; Veal *et al.*, 1988), nuclear

magnetic resonance (NMR) measurements (Viefhaus & Müller, 2006; Kim *et al.*, 2007) and single-crystal diffuse scattering analysis (Andersen *et al.*, 1986; Welberry *et al.*, 1992, 1995; Goff *et al.*, 1999). The latter has so far been restricted to selected layers in reciprocal space. More recent computational and technological advances allow the collection and interpretation of full three-dimensional data (Welberry *et al.*, 2003), opening up the possibility to re-visit controversially discussed defect clusters by utilizing the highly powerful 3D-ΔPDF technique (Weber & Simonov, 2012; Roth & Iversen, 2019; Simonov *et al.*, 2014), which represents differences in interatomic vector probabilities with respect to the average structure.

This work is organized as follows. We begin by presenting our experimentally obtained diffuse scattering data and quantitatively analyse the 3D-ΔPDFs for the dominant local interactions. We compare these results to literature findings and discuss the novel insights obtained by the 3D-ΔPDF analysis. Using our 3D-ΔPDF analysis we perform Monte Carlo simulations to show how the quantitative insight into local-order models can be realized in a model crystal. We conclude by discussing the significance of our quantitative 3D-ΔPDF analysis and outline the strengths of combining X-ray and neutron single-crystal diffuse measurements to solve complex local-order problems as encountered in YSZ.

## 2. Results and discussion

### 2.1. Sample material

The zirconia samples have a composition of $Zr_{0.82}Y_{0.18}O_{1.91}$, grown by the skull melting method, delivered by Djevahirdjan S. A., Monthey, Switzerland. The composition was confirmed by EDX measurements (see supporting information). For neutron measurements, the large clear single crystals were cut with a diamond saw into cubes with an edge length of approximately 5 mm. For X-ray diffraction measurements, the larger crystals were mechanically ground to a diameter of about 150 µm and polished. All measurements were carried out under ambient conditions.

### 2.2. X-ray and neutron diffraction experiments

X-ray diffraction experiments were performed on a Rigaku Synergy S diffractometer equipped with an Eiger 1M detector using Mo radiation. To avoid possible fluorescence a threshold of 17.4 keV was set for the detector. Simple $\varphi$ scans with 0.5° step width and 120 s exposure time were taken. 3D diffuse scattering data was reconstructed using the orientation matrix provided by CrysAlis PRO (Agilent, 2014) and custom Python scripts using Meerkat (Simonov, 2020). The experimentally obtained unit-cell parameter was $a = 5.1505 (5) \AA$.

Neutron diffraction experiments were carried out at the D19 instrument (λ = 0.95 Å, 0.1° steps, 80 s exposure per frame) (ILL, Grenoble, France) utilizing a 180° $\varphi$ scan. ILL are data available at https://dx.doi.org/10.5291/ILL-DATA.5-13-277. 3D diffuse scattering data reconstruction utilized the
orientation matrix as provided by *Int3d* (Katcho *et al.*, 2021) and a custom Python script.

The reflection conditions for space group *Fm$\overline{3}$m* were fulfilled in all cases and after careful inspection the data were symmetry averaged for *m$\overline{3}$m* Laue symmetry. Symmetry-averaged reciprocal space maps of selected layers are presented in Fig. 2.

## 2.3. 3D-ΔPDF maps

The general data processing procedure to obtain 3D-ΔPDF experiments is described by Koch *et al.* (2021). The experimentally obtained data were treated with the KAREN outlier rejection algorithm (Weng *et al.*, 2020). Additionally a custom punch-and-fill approach that interpolates the intensity in punched voxels was used to eliminate residual Bragg intensities. To avoid Fourier ripples the data were multiplied with a Gaussian falloff that smooths the edges of the measured reciprocal space section [see Weng *et al.* (2020)]. The FFT algorithm as implemented in *Meerkat* (Simonov, 2020) was used to obtain 3D-ΔPDF maps.

The 3D-ΔPDF maps in the *ab*0 layer and the *aac* layer are shown in Fig. 3. All 3D-ΔPDFs show a variety of different signatures. The most pronounced features are observed at the shortest interatomic vectors indicating strong local-order principles, which we will analyse in the following sections.

## 2.4. Oxygen nearest neighbour interactions

The (1/2, 0, 0) interatomic vectors are the shortest interatomic vectors in the idealized zirconia structure that only occur in the oxygen substructure. Detailed two-dimensional sections and three-dimensional renderings of the 3D-ΔPDFs around the (1/2, 0, 0) interatomic vector are shown in Fig. 4.

<figure><img src="image_7.png" /><figcaption>**Figure 2**
Reconstructed and symmetry-averaged diffuse scattering of YSZ obtained from neutron (left) and X-ray (right) diffraction experiments: (a) the *hk*0 layer and (b) the *hhl* layer.</figcaption></figure>

For the neutron diffraction experiments, a clear minimum shifted to shorter interatomic vectors by $\delta_{\overline{OO}}$ at (1/2 - $\delta_{\overline{OO}}$, 0, 0) and a clear maximum shifted to longer interatomic vectors by $\delta_{OO}^+$ at (1/2 + $\delta_{OO}^+$, 0, 0) can be observed. In the X-ray diffraction experiment the observed trend is the same but much weaker due to the much lower scattering power of oxygen as compared to the heavier metal ions. This is the typical signa-

<figure><img src="image_8.png" /><figcaption>**Figure 3**
3D-ΔPDF maps of YSZ obtained from neutron (left) and X-ray (right) diffraction experiments: (a) the *ab*0 layer and (b) the *aac* layer. Positive intensities in red, negative intensities in blue.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 4**
3D-ΔPDFs of YSZ obtained from neutron (left) and X-ray (right) diffraction experiments around the (1/2, 0, 0) interatomic vector. (a) Two-dimensional section in the *ab*0 layer. (b) Three-dimensional rendering of the intensity distribution. Volume shown in the region -0.15 ≤ *a* - 1/2, *b*, *c* ≤ 0.15. Black lines indicate the average structure interatomic vector at (1/2, 0, 0). 3D renderings are on the same relative intensity scale with respect to the minimum intensity in the rendering section. Positive intensities in red, negative intensities in blue.</figcaption></figure>
ture for a size-effect like relaxation (Weber & Simonov, 2012). A minimum in the 3D-ΔPDF indicates less scattering density than suggested by the average structure model, a maximum indicates more scattering density. The $\langle\frac{1}{2}, 0, 0\rangle$ interatomic vectors discussed here only occur in the oxygen substructure, hence only correlations that include oxygen ions and vacancies need to be considered. Any real structure configuration that includes a vacancy has zero scattering density, while any real structure configuration that involves two oxygen ions has more scattering density than suggested by the average structure ($b_O \times b_O > 0.955b_O \times 0.955b_O$). Hence we can attribute the minimum observed here to real structure configurations that involve a vacancy and occur with a higher probability at shorter interatomic vectors, while configurations that involve two oxygen ions can be attributed to the maximum and occur with a higher probability at longer interatomic vectors. This is a clear indication that locally oxygen ions relax towards neighbouring vacancies along the $\langle 1, 0, 0 \rangle$ directions.

This finding is consistent with computational and experimental reports from the literature, and qualitatively displayed in Fig. 1 (Frey *et al.*, 2005; Khan *et al.*, 1998; Fèvre *et al.*, 2005). The MD simulations (on $Zr_{0.9375}Y_{0.0675}O_{1.96875}$) by Fabris *et al.* (2002) report a relaxation of oxygen ions neighbouring a vacancy of 0.27 Å, the MD simulations (on $Zr_{0.865}Y_{0.135}O_{1.9325}$) of Fèvre *et al.* (2005) report a shift of 0.40 Å and the first-principle calculations (on $Zr_{0.9375}Y_{0.0675}O_{1.96875}$) of Stapper *et al.* (1999) report a shift of 0.24 Å. Experimental reports using Bragg data refinements from Goff *et al.* (1999) (on $Zr_{0.8}Y_{0.2}O_{1.9}$, shift 0.04 r.l.u. ≈ 0.20 Å), Kaiser-Bischoff *et al.* (2005) (on $Zr_{0.74}Y_{0.26}O_{1.87}$, shift 0.24 Å) and Ishizawa *et al.* (1999) (on $Zr_{0.758}Y_{0.242}O_{1.879}$, shift 0.31 Å) report similar results. These results show no clear indication of a correlation of the shift magnitude and the dopant concentration of the sample. Our results also support a clear shift of the oxygen ions along $\langle 1, 0, 0 \rangle$, and show no indication for previously suggested shifts along $\langle 1, 1, 1 \rangle$ directions (Ishizawa *et al.*, 1999; Argyriou *et al.*, 1996) as such a shift would distort the homogeneous maximum at $(\frac{1}{2} + \delta_{OO}^+, 0, 0)$. However, an additional, more isotropic off-axis relaxation is possible, as the observed maximum in the 3D-ΔPDF shows a disc-like feature. This indicates that for two oxygen ions that are separated by $(\frac{1}{2} + \delta_{OO}^a, \delta_{OO}^b, \delta_{OO}^c)$, $\delta_{OO}^a$ shows a much more narrow distribution than $\delta_{OO}^b$ and $\delta_{OO}^c$.

The 3D-ΔPDF analysis also allows a quantitative estimation of the shift magnitudes. For this purpose we fit the position of the minimum centred at $(\frac{1}{2} - \delta_{OO}^-, 0, 0)$ with a three-dimensional Gaussian distribution. The resulting parameter for $\delta_{OO}^-$ is $1.02 (1) \times 10^{-1}$ r.l.u. for the fit to the neutron data and $0.91 (1) \times 10^{-1}$ r.l.u. for the fit to the X-ray data. From electrostatic considerations we assume that neighbouring vacancies are highly unlikely. Therefore we interpret the minimum as generated by real structure oxygen-vacancy correlations. With the experimentally obtained unit-cell parameter $a = 5.1505 (5)$ Å and our refined shift magnitude we estimate that the oxygen ions neighbouring a vacancy relax 0.525 (5) Å along $\langle 1, 0, 0 \rangle$ towards the vacancy. The shift amplitude is larger than previous experimental reports that utilized Bragg

data analysis suggested and more similar to the MD simulations of Fèvre *et al.* (2005). Due to the much higher relative sensitivity of the neutron 3D-ΔPDF we used the parameter derived from the neutron data in the shift estimation. Nevertheless, the X-ray 3D-ΔPDF shows a clear signature at $\langle\frac{1}{2}, 0, 0\rangle$ and the quantitative analysis suggests that oxygen displacements can be directly observed in X-ray diffraction experiments, with a reasonable quantitative agreement to the neutron data.

## 2.5. Oxygen metal interactions

The $\langle\frac{1}{4}, \frac{1}{4}, \frac{1}{4}\rangle$ interatomic vectors only occur between the metal and the oxygen substructures. Detailed two- and three-dimensional 3D-ΔPDFs of the $(\frac{1}{4}, \frac{1}{4}, \frac{1}{4})$ interatomic vector are shown in Fig. 5 and resemble the typical signature for a positive atomic displacement parameter (ADP) correlation where the local bond distance variation is smaller than suggested by the ADPs of the average structure. The simple analysis and interpretation of interatomic vector as for the $\langle\frac{1}{2}, 0, 0\rangle$ vector is not possible here, as there are four distinct pair correlations that can be observed at the $\langle\frac{1}{4} \frac{1}{4} \frac{1}{4}\rangle$ interatomic vectors: possible pairs are Zr-O, Zr-vacancy, Y-O and Y-vacancy.

For X-ray diffraction, there is negligible contrast between $Zr^{4+}$ and $Y^{3+}$. For neutron diffraction the contrast between Zr and Y is not large either [$b(Zr) = 7.16$ versus $b(Y) = 7.75$], which does not allow for a direct disentanglement of signatures associated with Y-O and Zr-O distances. However, similar to the $\langle\frac{1}{2}, 0, 0\rangle$ interatomic vectors, any real space configuration that involves a vacancy has zero scattering density. Therefore, positive correlations indicate interatomic vectors with a higher probability for locally observed metal-oxygen vectors. In a similar fashion to the $\langle\frac{1}{2}, 0, 0\rangle$ interatomic

<figure><img src="image_10.png" /><figcaption>Figure 5</figcaption></figure>

3D-ΔPDFs obtained from neutron (left) and X-ray (right) diffraction experiments. (a) Two-dimensional section in the $ab0.25$ layer. (b) Three-dimensional rendering of the intensity distribution. Volume shown in the region $0.1 \le a, b, c \le 0.40$. Black lines indicate the average interatomic vector at $(\frac{1}{4}, \frac{1}{4}, \frac{1}{4})$. Positive intensities in red, negative intensities in blue.
vectors we fit the position of the maximum with a three-dimensional Gaussian distribution centred at $(\frac{1}{4} + \delta_{MO}, \frac{1}{4} + \delta_{MO}, \frac{1}{4} + \delta_{MO})$. The resulting parameters for $\delta_{MO}$ are $-1.58 (10) \times 10^{-2}$ r.l.u. for X-ray and $-0.95 (8) \times 10^{-2}$ r.l.u. for neutron diffraction. The observed maximum in the 3D-ΔPDF is shifted towards the centre of the 3D-ΔPDF space, contracted in the direction of the shift and elongated perpendicular of to the shift. This indicates a local contraction of the oxygen-metal bond length, which can be achieved if metal ions relax away from neighbouring vacancies along the (1, 1, 1) directions. The resulting local configurations show a variety of oxygen-metal bond distances, which accounts for the elongation of the maximum perpendicular to the shift direction. With the experimentally obtained unit-cell parameter $a = 5.1505 (5) \AA$, we estimate the average metal-O2- nearest neighbour distance at 2.09 (1) Å for X-ray and 2.15 (1) Å for neutron diffraction experiments. The differences in the calculated distances can be attributed to the differences in scattering contrasts. The fact that Y has the larger neutron scattering length and the estimated average metal-oxygen bond length is longer for the neutron refinement than for the X-ray refinement indicate that locally Y-O bond lengths are larger than Zr-O bond length, which is consistent with ionic radii (Shannon, 1976; Prince, 2004).

The arrangement of the oxygen ions and the vacancies in the oxygen substructure determines the coordination numbers of the metals. In the average structure of cubic ZrO₂, each Zr⁴⁺ ion is in regular cubic coordination, while monoclinic ZrO₂ only shows sevenfold coordination (Frey et al., 2005). Y₂O₃ shows sixfold coordinated metal ions at ambient conditions (Antic et al., 1993). Hence, it seems that both the Zr and the Y ions compete for lower coordination numbers; electrostatic effects and atomic sizes need to be taken into account when evaluating which type of metal ions neighbours a vacancy in the real structure. The simulations of Bogicevic & Wolverton (2003) and Khan et al. (1998) suggest that for YSZ NNN vacancies are preferred. Experimentally this is confirmed by the ⁸⁹Y solid-state MAS-NMR studies of Viefhaus & Müller (2006), which only show a significant onset of the peaks for six- and sevenfold coordinated Y ions above 10 mol% Y₂O₃, which is above the 9 mol% Y₂O₃ in our specimen.

Different metal coordination and different ionic radii lead to different oxygen-metal bond lengths. In Table 1 we summarize simulated bond distances given by Khan et al. (1998), experimentally measured bond distances from EXAFS studies in YSZ by Ishizawa et al. (1999) and Catlow et al. (1986) and the sum of effective ionic radii based on Shannon (1976). However, the presence of the oxygen vacancies and the associated relaxations are likely to yield not perfect coordination polyhedra but rather distorted bonding environments (see Fig. 1). Such distorted bond environments are also supported by density functional theory calculations on superstructures of Zr₁⁻ₓYₓO₂⁻ₓ/₂ (Pruned & Artacho, 2005), X-ray absorption studies (Li et al., 1993) and the fact that sevenfold coordination in pure ZrO₂ at room temperature shows a considerable variance in Zr⁴⁺-O²⁻ bond distances (between 2.05 Å and 2.26 Å) (Yashima et al., 1995).

<caption>Table 1</caption>

Metal-oxygen bond lengths (Å) as a function of coordination number (CN) and sum of effective ionic radii.

<table><thead><tr><th>Bond</th><th>CN</th><th>MD simulation†</th><th>EXAFS</th><th>Sum of effective ionic radii‡</th></tr></thead><tbody><tr><td rowspan="3">Y-O</td><td>6</td><td>2.267<sup>a</sup></td><td>2.32§</td><td>2.28</td></tr><tr><td>7</td><td>2.267<sup>a</sup></td><td>–</td><td>2.34</td></tr><tr><td>8</td><td>2.340<sup>b</sup></td><td>2.28¶</td><td>2.40</td></tr><tr><td rowspan="3">Zr-O</td><td>6</td><td>2.117<sup>b</sup></td><td>2.13§/2.11¶</td><td>2.10</td></tr><tr><td>7</td><td>2.117<sup>b</sup></td><td>2.13§/2.11¶</td><td>2.16</td></tr><tr><td>8</td><td>2.119<sup>a</sup></td><td>2.11¶</td><td>2.22</td></tr></tbody></table>

The estimated average metal-oxygen nearest neighbour distances from our 3D-ΔPDF analysis is, therefore, likely also an average obtained from several distorted bond environments. Comparing the deduced average distance with the sum of the atomic radii (Shannon, 1976) in Table 1 a sevenfold Zr coordination is supported. This is consistent with the EXAFS studies by Ishizawa et al. (1999) and Catlow et al. (1986), the ⁸⁹Y solid-state MAS-NMR studies of Viefhaus & Müller (2006), as well as the computational work of Bogicevic & Wolverton (2003) and Khan et al. (1998) that favour dopant ions as NNN to the vacancies, which results in a preference for Zr in sevenfold coordination.

## 2.6. Metal-metal interactions

The ⟨1/2, 1/2, 0⟩ are the most complex of the nearest neighbour interatomic vectors as they occur in both structures. The maxima associated with O-O correlations are overlaid by metal-metal interactions. 3D-ΔPDFs of the (1/2, 1/2, 0) interatomic vector are shown in Fig. 6.

<figure><img src="image_11.png" /><figcaption>Figure 6</figcaption></figure>

3D-ΔPDFs obtained from neutron (left) and X-ray (right) diffraction experiments. (a) Two-dimensional plot in the ab0 layer. (b) Three-dimensional rendering of the intensity distribution. Volume shown in the region 0.25 ≤ a, b, c + 0.5 ≤ 0.75. Black lines indicate the average interatomic vector at (1/2, 1/2, 0). Positive intensities in red, negative intensities in blue.

<footnote>† Khan et al. (1998). The MD simulation distinguishes between oxygen vacancy NN to dopant ion and oxygen vacancy NNN to dopant ion. We interpret this in terms of coordination numbers. (a) NN dopant vacancy, (b) NNN dopant vacancy. § Ishizawa et al. (1999). ¶ Catlow et al. (1986). ‡ Shannon (1976).</footnote>
For both our diffraction experiments, the observed signatures are complex and the differences highlight the contrast between X-ray and neutron diffraction experiments. However, both signatures show two distinct locations of maxima, which can be interpreted as a result of different local configurations: configurations where both or one of the bridging oxygen ions are missing in between the metals will result in a different metal-metal interatomic vector than those of configurations where both bridging oxygen ions are present. Furthermore the exact intensity and distance distribution will depend on the coordination numbers and types of metals involved. The entanglement of the interatomic distances involved in the signature at the (1/2, 1/2, 0) vectors is therefore effectively a multi body correlation and can only be accessed indirectly with scattering experiments where only pair correlations can be probed directly (Welberry, 2004; Baake & Grimm, 2009). This limitation is also present in 3D-ΔPDF and a direct interpretation of the signature at (1/2, 1/2, 0) relies on information deduced from shorter interatomic vectors.

However, assuming a shift of the metal ions neighbouring a vacancy away from the vacancy, as deduced from our data in the previous section, we can assign the maximum observed at (1/2 + δ<sub>MM</sub>, 1/2 + δ<sub>MM</sub>, 0): for simplicity we consider a single pair of metal ions M1 at (0, 0, 0) and M2 at (1/2, 1/2, 0) with the bridging vacancy located at (1/4, 1/4, 1/4). The metal ions will relax away from the vacancies in the (1, 1, 1) directions. M1 will be shifted to (−Δ<sub>MO</sub>, −Δ<sub>MO</sub>, −Δ<sub>MO</sub>) and M2 will be shifted to (1/2 + Δ<sub>MO</sub>, 1/2 + Δ<sub>MO</sub>, −Δ<sub>MO</sub>), yielding the new inter-atomic vector between M1 and M2 as (1/2 + 2Δ<sub>MO</sub>, 1/2 + 2Δ<sub>MO</sub>, 0). While the shift δ<sub>MO</sub> determined in the previous section yields an average metal-oxygen bond length, the shift Δ<sub>MO</sub> we determine here is directly related to the shift of a metal atom neighbouring a vacancy. We quantify this shift by fitting the maximum at (1/2 + δ<sub>MM</sub>, 1/2 + δ<sub>MM</sub>, 0) in the X-ray 3D-ΔPDF

with a three-dimensional Gaussian distribution. We choose the X-ray data for this estimation as the scattering contributions here are dominated by the contributions of the metal ions. The fitted shift magnitude is δ<sub>MM</sub> = 2Δ<sub>MO</sub> = 1.043 (4) × 10<sup>-1</sup> r.l.u. This corresponds to a metal shift along the (1, 1, 1) directions away from the vacancy of 0.465 (2) Å, which is slightly more pronounced than references in the literature. The MD simulations of Fabris *et al.* (2002) (Zr<sub>0.9375</sub>Y<sub>0.0675</sub>O<sub>1.96875</sub>) report a metal shift away from the neighbouring vacancy of 0.18 Å, the MD simulations of Fèvre *et al.* (2005) (Zr<sub>0.865</sub>Y<sub>0.135</sub>O<sub>1.9325</sub>) report a shift of 0.1 Å and the first-principle calculations of Stapper *et al.* (1999) (Zr<sub>0.9375</sub>Y<sub>0.0675</sub>O<sub>1.96875</sub>) report a shift of 0.18 Å. Experimental reports using Bragg data refinements from Goff *et al.* (1999) (Zr<sub>0.8</sub>Y<sub>0.2</sub>O<sub>1.9</sub>, shift of 0.028 r.l.u. ≈ 0.11 Å), Kaiser-Bischoff *et al.* (2005) (Zr<sub>0.74</sub>Y<sub>0.26</sub>O<sub>1.87</sub>, shift 0.22 Å) and Ishizawa *et al.* (1999) (Zr<sub>0.758</sub>Y<sub>0.242</sub>O<sub>1.879</sub>, 0.0219 r.l.u. shift 0.08 Å) vary largely in the reported shift magnitude and show no clear indication of a correlation of the shift magnitude and the dopant concentration of the sample.

## 2.7. Possibility of vacancy clustering and interstitial metal atoms at (1/2, 1/2, 1/2)

Questions that have been discussed controversially in the literature are the possibility of vacancy clusters along the (1/2, 1/2, 1/2) vectors – eventually forming chains of pyrochlore-like structural elements (Welberry *et al.*, 1993; Goff *et al.*, 1999) – and the possibility of interstitial metal ions at (1/2, 1/2, 1/2) (Goff *et al.*, 1999). Both of these questions can be clarified by analysing the 3D-ΔPDFs in the vicinity of (1/2, 1/2, 1/2) (see Fig. 7).

The possibility of metal ion interstitials at (1/2, 1/2, 1/2) as considered by Goff *et al.* (1999) can also be ruled out here by the absence of relevant signatures in the X-ray 3D-ΔPDF in Fig. 7. Possible interstitial metal ions at (1/2, 1/2, 1/2) would lead to metal-metal interatomic vectors at (1/2, 1/2, 1/2), which are not described by our average structural model and would hence result in a signature in the 3D-ΔPDF. The lower scattering power of oxygen compared to the metal ions – supported by the relative strength of the signatures in the X-ray 3D-ΔPDFs observed at (1/2, 0, 0) and (1/2, 1/2, 0) – suggests that the X-ray 3D-ΔPDF would be very sensitive to such a signature, while in the neutron case the signature present at (1/2, 1/2, 1/2) can be interpreted in terms of oxygen-oxygen correlations.

The shortest oxygen-oxygen interatomic vectors are (1/2, 0, 0) and (1/2, 1/2, 0). For these vectors vacancy clustering becomes highly unlikely from an electrostatic viewpoint (Bogicevic & Wolverton, 2003). The shortest vectors where vacancy pairs are allowed from electrostatic considerations are the (1/2, 1/2, 1/2) vectors. This yields the possibility to obtain six-coordinated metal ions in the structure. The doping level in the specimen at hand (Zr<sub>0.82</sub>Y<sub>0.18</sub>O<sub>1.91</sub>) is low enough to distribute the vacancies in the structure without needing to form sixfold coordinated ions, unlike for higher-doped variants such as described by Welberry *et al.* (1992). The in-depth simulation-based study of Bogicevic & Wolverton

<figure><img src="image_12.png" /><figcaption>**Figure 7** 3D-ΔPDFs obtained from neutron (left) and X-ray (right) diffraction experiments. (a) Two-dimensional plot in the ab0.5 layer. (b) Three-dimensional rendering of the intensity distribution. Volume shown in the region 0.25 ≤ a, b, c ≤ 0.75. Black lines indicate the average interatomic vector at (1/2, 1/2, 1/2). Positive intensities in red, negative intensities in blue.</figcaption></figure>
(2003), however, suggests that the formation of double vacancies along $\langle\frac{1}{2}, \frac{1}{2}, \frac{1}{2}\rangle$ can nevertheless be favourable.

The signatures in Fig. 7 are very weak in the case of X-ray diffraction – basically indistinguishable from residual background noise. The situation is different for the neutron scattering case: here the scattering length of oxygen is comparable to that of the metals and the relative strength of signatures attributed to oxygen-oxygen interactions is much greater than in the case of X-ray diffraction. Here, we observe a clear positive signature at $(\frac{1}{2}, \frac{1}{2}, \frac{1}{2})$ suggesting a positive correlation. If we consider pure substitutional disorder, positive correlations indicate a higher probability for likewise neighbours at this interatomic vector (Weber & Simonov, 2012). In our case this indicates a higher probability for oxygen-oxygen and vacancy-vacancy pairs along $\langle\frac{1}{2}, \frac{1}{2}, \frac{1}{2}\rangle$, i.e. a tendency to form double-vacancy pairs along this interatomic vector. Comparing the strength of the signature with the signature at $\langle\frac{1}{2}, 0, 0\rangle$ (note Figs. 4 and 7 are on the same relative scale), the signature of the double vacancies along $\langle\frac{1}{2}, \frac{1}{2}, \frac{1}{2}\rangle$ is much weaker. This suggests that there is no strict formation of double vacancies which would in turn lead to formation of zigzag chains of six-coordinated metal ions and local pyrochlore-like structural elements. This observation is in agreement with the findings of Welberry et al. (1995).

## 2.8. Longer-range interactions and extent of correlations

A direct interpretation of longer-range interactions in terms of underlying configurations is basically impossible due to the configurational complexity of possible arrangements of nearest neighbour, next nearest neighbour and next next nearest neighbour oxygen-vacancy pairs and metal-metal pairs. Nevertheless, the obtained 3D-ΔPDFs in Fig. 3 allow a direct judgement of the extent to which correlations due to the previously identified local correlations persist. For both the neutron and X-ray 3D-ΔPDF we do not observe correlations that are significantly longer than two unit cells, i.e. ≈ 10 Å, which is well below the estimated observable limit (see supporting information) and significantly shorter than previously reported in literature (Goff et al., 1999; Frey et al., 2005).

Possible differences in the neutron and X-ray 3D-ΔPDFs allow distinguishing whether the correlations are propagated in both substructures to the same extent: the limited sensitivity of X-ray 3D-ΔPDFs to oxygen-oxygen correlations allows the interpretation of the extent of correlations in terms of metal displacements only. In the X-ray 3D-ΔPDF, correlations along $\langle 1, 1, 0 \rangle$ persist longer than in the neutron case (see e.g. signature at $\langle \frac{3}{2}, \frac{3}{2}, 0 \rangle$ in Fig. 3), providing evidence that the metal-metal nearest neighbour distortion described in §2.6 is propagated to further neighbours in the same direction.

## 2.9. 3D-ΔPDF informed modelling

The detailed analysis of the experimentally obtained 3D-ΔPDFs we performed in the previous sections lays the groundwork for the simulation of a simplistic atomistic model that realizes the correlations we derived. For this purpose we

use three successive Monte Carlo simulations, where the first two simulations establish the chemical ordering and the third simulation relaxes the atomistic positions according to our analysis. We simulate five model crystals of $10 \times 10 \times 10$ unit cells. First the metal ions are distributed at random to match the average composition of 82% Zr and 18% Y ions. Charge balance is ensured by removing the suitable amount of oxygen ions.

The first Monte Carlo simulations induces chemical ordering of the metal ions. We assume for electrostatic reasons that Y ions tend to avoid being close to each other. The Monte Carlo simulations swap Zr and Y ions to avoid NN and NNN Y-Y pairs.

<figure><img src="image_13.png" /><figcaption>Figure 8</figcaption></figure>

3D-ΔPDFs obtained from our simple atomistic model (M) from neutron (left) and X-ray (right) diffuse scattering calculations compared to the experimentally obtained 3D-ΔPDFs (E). (a) Two-dimensional section around $(\frac{1}{2}, 0, 0)$ in the ab0 layer. Model top, experiment bottom. (b) Two-dimensional section around $(\frac{1}{4}, \frac{1}{4}, \frac{1}{2})$ in the ab0.25 layer. Model top left, experiment bottom right. (c) Two-dimensional section around $(\frac{1}{2}, \frac{1}{2}, 0)$ in the ab0 layer. Model top left, experiment bottom right. (d) Two-dimensional section around $(\frac{1}{2}, \frac{1}{2}, \frac{1}{2})$ in the ab0.5 layer. Model top left, experiment bottom right.
The second Monte Carlo simulation then swaps oxygen ions and vacancies in such a fashion that:

(i) Y-ions prefer eightfold coordination and avoid being next to vacancies,

(ii) Vacancy pairs separated by $\langle\frac{1}{2}, 0, 0\rangle$ and $\langle\frac{1}{2}, \frac{1}{2}, 0\rangle$ are penalized with a high energy,

(iii) A minimal energy gain is introduced for vacancies separated by $\langle\frac{1}{2}, \frac{1}{2}, \frac{1}{2}\rangle$.

The resulting structures fulfil (i) and (ii) without violations. On average 3.06 (13)% of the metal ions are sixfold coordinated.

The third Monte Carlo simulation then introduces relaxations as deduced from the analysis of the 3D-ΔPDFs: assuming that every vacancy introduced in the structure causes its six neighbouring oxygen ions to relax towards the vacancy along the $\langle 1, 0, 0 \rangle$ directions and its four nearest neighbour metal ions to relax along $\langle 1, 1, 1 \rangle$ directions away from the vacancy, we introduced a static shift for the respective percentage of the ions and then used the Monte Carlo algorithm as implemented in the *DISCUS* program (Neder & Proffen, 2008) to switch the displacements. The energy targets were set as spring potentials with metal-oxygen target distances taken as the sum of the ionic radii (Shannon, 1976).

The diffuse scattering (see supporting information) was calculated using the *DISCUS* program (Neder & Proffen, 2008) on a grid adapted to the size of the supercell. To obtain the model 3D-ΔPDFs we used a customized punch and fill algorithm (punch size was one voxel, as there is no experimental broadening of Bragg reflections in the calculated diffuse scattering). The same Gaussian falloff was multiplied to the simulated data as was to the experimental data and *Meerkat* (Simonov, 2020) was used to obtain the resulting simulated 3D-ΔPDFs in the regions of interest shown in Fig. 8 and in the supporting information.

A comparison of the 3D-ΔPDFs from our simple model and the respective sections in the experimentally obtained 3D-ΔPDFs in Fig. 8 shows that the nature of the major features that we quantified in our analysis are reproduced. The majority of the positions of maxima and minima are reproduced well while we observe more significant differences in the intensity distribution and the spatial extent of the features. These discrepancies we attribute to the simplicity of our model that only contains the quantities derived from a direct interpretation of the experiment. A refinement, e.g. using a reverse Monte Carlo simulation, would likely result in a much better agreement between data and model. However, our emphasis here is on the direct interpretation of the 3D-ΔPDF and our model confirms that the quantities we derive from this direct interpretation can be used to build a simple model that reproduces the nature of the local interactions and therefore confirms the results from our direct interpretation.

### 3. Concluding remarks

The defect structure of YSZ has been previously studied by several computational and experimental techniques (Frey *et al.*, 2005). Bragg data analysis, NMR and EXAFS studies

provide only limited information about the local order in the system and a comprehensive model of the defect structure can only be obtained from single crystal diffuse scattering analysis. In the past, such analysis required computationally expensive modelling, and a direct interpretation of present distortions from the measured reciprocal space sections required very extensive expert knowledge or was simply impossible (Andersen *et al.*, 1986; Welberry *et al.*, 1993, 1995; Goff *et al.*, 1999). With the advances in experimental techniques it is now possible to obtain full three-dimensional reciprocal space coverage and the application of the 3D-ΔPDF allows direct interpretation of the data in terms of defect models (Weber & Simonov, 2012; Roth & Iversen, 2019; Simonov *et al.*, 2014).

In our contribution, we demonstrated how to combine the information from X-ray and neutron 3D-ΔPDF to directly establish a quantitative defect model for YSZ. The signatures we observe on the shortest interatomic vectors, i.e. $\langle\frac{1}{2}, 0, 0\rangle$, $\langle\frac{1}{4}, \frac{1}{4}, \frac{1}{4}\rangle$, $\langle\frac{1}{2}, \frac{1}{2}, 0\rangle$ and $\langle\frac{1}{2}, \frac{1}{2}, \frac{1}{2}\rangle$, enabled us to build a conclusive model for local correlations. The combination of the two different types of radiation is crucial for the quantitative analysis: the X-ray 3D-ΔPDF has a limited sensitivity to oxygen-oxygen correlations, while in the neutron 3D-ΔPDF the signatures of metal-metal correlations and oxygen-oxygen correlations overlap. The combined analysis allows the signatures to be disentangled and permits the quantitative analysis of local relaxation mechanisms. The shift directions we proposed from our analysis are in good agreement with previously reported relaxations from various experimental investigations and simulation efforts, the shift magnitudes we analyse are more pronounced than previously suggested. We consider our work as an important first step in the progress towards a direct methods type analysis for diffuse scattering using 3D-ΔPDF methods.

Our study shows how a simple, direct interpretation of the 3D-ΔPDF leads to results that align well with established interpretations in literature, that involved complex modelling of the described correlations. In agreement with (Goff *et al.*, 1999) we find no evidence for interstitial metal atoms at $\langle\frac{1}{2}, \frac{1}{2}, \frac{1}{2}\rangle$. We find strong evidence for double vacancy pairs along $\langle\frac{1}{2}, \frac{1}{2}, \frac{1}{2}\rangle$ (Bogicevic & Wolverton, 2003; Welberry *et al.*, 1993). Our analysis confirms oxygen displacements along $\langle 1, 0, 0 \rangle$ and metal displacements along $\langle 1, 1, 1 \rangle$ (Kaiser-Bischoff *et al.*, 2005), while we do not find direct evidence for oxygen displacements along $\langle 1, 1, 1 \rangle$ (Ishizawa *et al.*, 1999; Argyriou *et al.*, 1996) or metal displacements along $\langle 1, 1, 0 \rangle$ (Welberry *et al.*, 1993).

We present a simplistic static local-order model that describes the most local interatomic correlations well, but does not reproduce all the signatures observed in the data; in particular, the observed longer-range correlations are not described in the model presented here. One way of incorporating such longer-range interactions would be a big box modelling approach, e.g. a reverse Monte Carlo simulation: by fixing the directions and the magnitudes of the shifts to the values identified in our 3D-ΔPDF analysis, a reverse Monte Carlo simulation would be able to generate a model crystal that captures longer-range correlations that are not directly
interpretable from the 3D-ΔPDF analysis. This in turn would provide a tool to identify longer-range interactions that are not present in the current model.

Our measurements are performed at ambient conditions; hence the model only captures correlations that are present at ambient conditions. However, for the technical application of YSZ as an oxygen ion conductor, the local order at elevated temperatures may be more relevant (Devanathan et al., 2006; Kaiser-Bischoff et al., 2005; Krishnamurthy et al., 2004; Itoh et al., 2015). For variable-temperature studies, the 3D-ΔPDF analysis as we present here facilitates a direct and quantitative comparison of correlations in real space as a function of temperature, providing a more complete picture than previously analysed selected sections of reciprocal space.

Another aspect that is of importance for technological applications is the interplay of compositional disorder and lattice dynamics. The picture presented here only describes a static local-order model and therefore does not cover the full picture of correlations at ambient conditions. The possible coupling of the local order to the lattice dynamics of a system is a particular challenge (Ziman, 1979; Snyder & Toberer, 2011). In YSZ the phonon anharmonicity at lower temperatures has been attributed to the defect structure and in turn to a higher vibrational entropy (Li et al., 2015). A comprehensive, quantitative and reliable local-order model is needed for a complete understanding of the lattice dynamics in complex systems such as YSZ. In the future, the combination of temperature dependent energy discriminated diffraction experiments and supercell lattice dynamical calculations (Overy et al., 2017; Schmidt et al., 2022) based on our local-order model will enable a complete understanding of the structure and dynamics of such a complex locally ordered system.

## Acknowledgements

We are grateful to Colin Johnston (Oxford University) for running and interpreting EDX measurements. Nebil Ayape Katcho (ILL) is thanked for the detailed help with the neutron data reconstruction. Arkadiy Simonov (ETH Zürich) is acknowledged for providing the *Meerkat* software. Yahua Liu, Feng Ye, Christina Hoffmann (all Oak Ridge National Laboratory) and Paul Benjamin Klar (University of Bremen) are acknowledged for fruitful discussions. Financial support for James D. Martin through the Gastprofessorenprogramm at the Friedrich-Alexander-Universität is gratefully acknowledged. Open access funding enabled and organized by Projekt DEAL.

## Funding information

Funding for this research was provided by: H2020 European Research Council (grant No. 788144 to Andrew L Goodwin).

## References

* Agilent (2014). *CrysAlis PRO*. Agilent Technologies Ltd, Yarnton, Oxfordshire, England.

* Andersen, N. H., Clausen, K., Hackett, M., Hayes, W., Hutchings, M., Macdonald, J. E. & Osborn, R. (1986). *Phys. B+C*, **136**, 315–317.

* Antic, B., Önnerud, P., Rodic, D. & Tellgren, R. (1993). *Powder Diffr.* **8**, 216–220.

* Argyriou, D. N., Elcombe, M. M. & Larson, A. C. (1996). *J. Phys. Chem. Solids*, **57**, 183–193.

* Baake, M. & Grimm, U. (2009). *Phys. Rev. B*, **79**, 020203.

* Bogicevic, A. & Wolverton, C. (2003). *Phys. Rev. B*, **67**, 024106.

* Bondars, B., Heidemane, G., Grabis, J., Laschke, K., Boysen, H., Schneider, J. & Frey, F. (1995). *J. Mater. Sci.* **30**, 1621–1625.

* Boysen, H., Frey, F. & Vogt, T. (1991). *Acta Cryst.* **B47**, 881–886.

* Catlow, C., Chadwick, A., Greaves, G. & Moroney, L. (1986). *J. Am. Ceram. Soc.* **69**, 272–277.

* Devanathan, R., Weber, W. J., Singhal, S. C. & Gale, J. D. (2006). *Solid State Ionics*, **177**, 1251–1258.

* Fabris, S. (2002). *Acta Mater.* **50**, 5171–5178.

* Fèvre, M., Finel, A. & Caudron, R. (2005). *Phys. Rev. B*, **72**, 104117.

* Frey, F., Boysen, H. & Kaiser-Bischoff, I. (2005). *Z. Kristallogr. Cryst. Mater.* **220**, 1017–1026.

* Goff, J., Hayes, W., Hull, S., Hutchings, M. & Clausen, K. N. (1999). *Phys. Rev. B*, **59**, 14202–14219.

* Hund, F. (1951). *Z. Elektrochem. Angew. Phys. Chem.* **55**, 363–366.

* Ishizawa, N., Matsushima, Y., Hayashi, M. & Ueki, M. (1999). *Acta Cryst.* **B55**, 726–735.

* Itoh, T., Mori, M., Inukai, M., Nitani, H., Yamamoto, T., Miyanaga, T., Igawa, N., Kitamura, N., Ishida, N. & Idemoto, Y. (2015). *J. Phys. Chem. C*, **119**, 8447–8458.

* Kaiser-Bischoff, I., Boysen, H., Scherf, C. & Hansen, T. (2005). *Phys. Chem. Chem. Phys.* **7**, 2061–2067.

* Katcho, N. A., Cañadillas-Delgado, L., Fabelo, O., Fernández-Díaz, M. T. & Rodríguez-Carvajal, J. (2021). *Crystals*, **11**, 897.

* Khan, M. S., Islam, S. M. & Bates, D. R. (1998). *J. Mater. Chem.* **8**, 2299–2307.

* Kim, N., Hsieh, C.-H., Huang, H., Prinz, F. B. & Stebbins, J. F. (2007). *Solid State Ionics*, **178**, 1499–1506.

* Koch, R. J., Roth, N., Liu, Y., Ivashko, O., Dippel, A.-C., Petrovic, C., Iversen, B. B., Zimmermann, M. v., & Bozin, E. S. (2021). *Acta Cryst.* **A77**, 611–636.

* Krishnamurthy, R., Yoon, Y.-G., Srolovitz, D. & Car, R. (2004). *J. Am. Ceram. Soc.* **87**, 1821–1830.

* Li, C. W., Smith, H. L., Lan, T., Niedziela, J. L., Muñoz, J. A., Keith, J. B., Mauger, L., Abernathy, D. L. & Fultz, B. (2015). *Phys. Rev. B*, **91**, 144302.

* Li, P., Chen, I.-W. & Penner-Hahn, J. E. (1993). *Phys. Rev. B*, **48**, 10074–10081.

* Morinaga, M., Cohen, J. B. & Faber, J. (1979). *Acta Cryst.* **A35**, 789–795.

* Neder, R. B. & Proffen, T. (2008). *Diffuse Scattering and Defect Structure Simulations: a cook book using the program DISCUS*, IUCr Texts on Crystallography, No. 11. Oxford University Press.

* Overy, A. R., Simonov, A., Chater, P. A., Tucker, M. G. & Goodwin, A. L. (2017). *Phys. Status Solidi B*, **254**, 1600586.

* Pham, A. Q. & Glass, R. S. (1998). *Electrochim. Acta*, **43**, 2699–2708.

* Prince, E. (2004). Editor. *International Tables for Crystallography, Vol. C, Mathematical, Physical and Chemical Tables*, 3rd ed. Dordrecht: Springer.

* Pruneda, J. & Artacho, E. (2005). *Phys. Rev. B*, **72**, 085107.

* Roth, N. & Iversen, B. B. (2019). *Acta Cryst.* **A75**, 465–473.

* Schindler, K., Schmeisser, D., Vohrer, U., Wiemhöfer, H. & Göpel, W. (1989). *Sens. Actuators*, **17**, 555–568.

* Schmidt, E. M., Thomas, S., Bulled, J. M., Minelli, A. & Goodwin, A. L. (2022). *Acta Cryst.* **B78**, 385–391.

* Shannon, R. D. (1976). *Acta Cryst.* **A32**, 751–767.

* Simonov, A. (2020). *Meerkat*. Version 0.3.7. A program for reciprocal space reconstruction. https://github.com/aglie/meerkat.

* Simonov, A., Weber, T. & Steurer, W. (2014). *J. Appl. Cryst.* **47**, 2011–2018.
Snyder, G. J. & Toberer, E. S. (2011). In *Materials for sustainable energy: a collection of peer-reviewed research and review articles from Nature Publishing Group*, pp. 101–110. World Scientific.

Stapper, G., Bernasconi, M., Nicoloso, N. & Parrinello, M. (1999). *Phys. Rev. B*, **59**, 797–810.

Tsampas, M., Sapountzi, F. & Vernoux, P. (2015). *Catal. Sci. Technol.* **5**, 4884–4900.

Veal, B., McKale, A., Paulikas, A., Rothman, S. & Nowicki, L. (1988). *Phys. B+C*, **150**, 234–240.

Viefhaus, T. & Müller, K. (2006). *Int. J. Mater. Res.* **97**, 508–512.

Weber, T. & Simonov, A. (2012). *Z. Kristallogr.* **227**, 238–247.

Welberry, T., Butler, B., Thompson, J. & Withers, R. (1993). *J. Solid State Chem.* **106**, 461–475.

Welberry, T. R., Goossens, D. J., Haeffner, D. R., Lee, P. L. & Almer, J. (2003). *J. Synchrotron Rad.* **10**, 284–286.

Welberry, T., Withers, R. & Mayo, S. (1995). *J. Solid State Chem.* **115**, 43–54.

Welberry, T., Withers, R., Thompson, J. & Butler, B. (1992). *J. Solid State Chem.* **100**, 71–89.

Welberry, T. R. (2004). *Diffuse X-ray Scattering and Models of Disorder. IUCr Monographs on Crystallography*, Vol. 16. Oxford University Press.

Weng, J., Dill, E. D., Martin, J. D., Whitfield, R., Hoffmann, C. & Ye, F. (2020). *J. Appl. Cryst.* **53**, 159–169.

Yashima, M., Hirose, T., Katano, S., Suzuki, Y., Kakihana, M. & Yoshimura, M. (1995). *Phys. Rev. B*, **51**, 8018–8025.

Ziman, J. M. (1979). *Models of Disorder. The Theoretical Physics of Homogeneously Disordered Systems*. Cambridge University Press.