Tuneable local order in thermoelectric crystals

Distinct ‘hidden’ phases of a technologically relevant thermoelectric material, which are identical in terms of composition and periodic crystal structure, but differ on a local scale, are observed and can be controlled through synthesis conditions. The local structure is explained in terms of a vacancy repulsion model, and relaxations around vacancies are characterized.

behaviour. The average periodic structure gives rise to sharp peaks in a diffraction experiment (Bragg peaks), which can be readily analysed. The local deviations give rise to weak diffuse scattering, which is much more difficult to measure and especially to interpret on a structural basis (1-3).
Defective half-Heusler materials X1-xYZ, such as Nb1-xCoSb, have excellent thermoelectric properties and they show signs of local atomic order different from their average periodic order (4)(5)(6)(7)(8). Electron diffraction data have revealed structured diffuse scattering, which differ from sample to sample (9). The differences were first explained as a result of different nominal sample stoichiometry (9), but recently the differences in shape of the reported diffuse scattering has been theoretically explained by a simple model for effective vacancy repulsion on the disordered X substructure, without invoking different sample stoichiometry (10). The model suggests that the structure places vacancies as far apart as possible, essentially giving the effect of vacancy repulsion. Importantly, the model predicts that the degree of local order can be influenced by synthesis conditions, such as thermal quenching. Experimental tuning of local order would provide a new handle in materials research, which may allow for controlling properties in disordered systems.
Recently, Goodwin and coworkers used diffuse x-ray scattering data to establish vacancy distributions in Prussian Blue Analogue (PBA) materials (3). Different PBAs crystals showed dissimilar diffuse scattering patterns, and indeed PBA materials are known to have large variation e.g. in battery properties (11).
However, an understanding of how local structure can be controlled and how the local structure correlates with material properties is still largely unexplored. Here, we show that the synthesis conditions change the degree of local order in thermoelectric Nb1-xCoSb. Using diffuse synchrotron x-ray scattering data measured on single crystals, we first validate the theoretical vacancy repulsion model, and then we model the structural relaxation around the vacancies to provide direct experimental quantification of the local structure in these systems.

Results
The theoretical model for vacancy repulsion in defective half-Heusler materials predicts that samples quenched (Q) from high temperature will have a different degree of local order than slowly cooled (SC) samples (10). Using induction furnaces, two distinct types of samples have been synthesized, quenched samples from a levitation-melt procedure and slowly cooled samples using crucibles. Samples were prepared using both methods with nominal stoichiometries of Nb0.81CoSb and Nb0.84CoSb, and they are named such that "Q-0.81" is a quenched sample with nominal stoichiometry Nb0.81CoSb.
From very high-quality single-crystal x-ray scattering measurements performed at a synchrotron we obtain the average crystal structure for each single crystal from the Bragg peaks, including an accurate stoichiometry, In addition, from the same measurement we obtain the diffuse scattering, allowing for analysis of the local order. With this method we can directly correlate the diffuse scattering and local order with the periodic crystal structure and sample composition, as this is all obtained from one measurement for each single crystal. It further avoids errors in determining the sample composition using other methods, as impurity phases are known to occur in these compounds (9).

Average Crystal structure
The average periodic crystal structures in all samples were quantified by analysis of the Bragg diffraction peaks. All samples have cubic average structures in space group 4 � 3 with the cell lengths in a narrow range between 5.894 and 5.899 Å at 300 K. The quality of the Bragg data is excellent with internal agreement (Rint) between 4.0% to 5.4% even for data resolutions with dmin < 0.43 Å, and average redundancies of 19-36 (see Table 1). Previous studies describe the average structure as an ideal half-Heusler structure with vacancies on the Nb sites (7). The ideal half-Heusler structure (space group 4 � 3 ) has Nb1-x at [0, 0, 0], Sb at [1/2, 1/2, 1/2] and Co at [1/4, 1/4, 3/4]. Nb1-x and Sb form a rock-salt structure, while Nb and Co form the sphalerite structure, as illustrated in Fig. 1a. Refining a model with free Nb occupancy gives good (low) agreement factors (R1: 2.29% -3.23%, wR2: 5.94% -7.39%), which normally would be considered excellent fit, especially with the large range of data. However, clear residuals are observed around the Sb and Co sites suggesting disorder, as shown in detail in the supplementary materials. If Sb is positioned off-center at [½+Δ, ½, ½] with refinement of Δ, and Co is offcentered at [1/4-δ, 1/4-δ, 3/4-δ] with refinement of δ, the agreement factors improve significantly for all samples (R1: 0.49% -1.06%, wR2: 0.78% -1.33%), Table 1. In this model, each Sb is split into six positions (each with 1/6 occupancy) forming an octahedron with corners pointing towards the six neighbouring Nb1-x sites, while Co is split into a tetrahedron with corners pointing towards neighbouring Sb, see Fig. 1b. For all samples, the Sb shift refines to Δ ~ 0.025 corresponding to a movement of 0.144 Å, and Co shifts δ ~0.0127 corresponding to a movement of 0.130 Å. Interestingly, the refined Nb occupancies are in the narrow range 0.820 to 0.827 for the five samples, suggesting that all Nb1-xCoSb samples are very close in stoichiometry with x ~ 1/6. Thus, there is no correlation between the nominal and refined sample stoichiometry. It is remarkable the average structure of all samples is identical with x ~ 1/6 irrespective of synthesis method and nominal sample stoichiometry. In all samples, Sb is shifted off-center by around 0.144 Å and Co by around 0.130 Å.
From refinements of the Bragg data all samples are identical, but "hidden" local structures are exposed through analysis of the diffuse scattering.  The measured x-ray scattering data, Fig. 2a, differ significantly from the electron diffraction data reported by Xia et al. (9), where the diffuse scattering consists of rings in the H0L plane with approximately constant intensity around the perimeter. In the x-ray data, the rings have clear intensity modulations. The difference is due to strong multiple scattering in the electron diffraction data, where the electron beam of strong Bragg peaks is re-scattered to give diffuse scattering averaged over Brillouin zones. This creates rings of approximately constant intensity, especially when data is measured along the zone-axis where Bragg scattering is strong as done by Xia et al. (9) The theoretical model for the vacancy distribution (10) was inspired by the electron diffraction data, and it produces calculated rings without the intensity modulation observed in the present x-ray scattering data. As will be shown, the intensity modulation of the diffuse scattering is mainly the result of structural relaxation of Sb and Co around vacancies, but overall the vacancy ordering model is validated. Similar modulations due to structural relaxation around vacancies, a type of size-effect, have been reported in other compounds (12,13).
To get a direct view of the local correlations in the samples, the diffuse scattering is Fourier transformed to obtain the thee-dimensional difference pair distribution function (3D-ΔPDF), which is the autocorrelation of the difference electron density: is the difference between the total electron density of the crystal and the periodic average electron density. The 3D-ΔPDF gives a direct view of the local deviations from the periodic average structure (14)(15)(16). Positive/negative features show vectors for which the real structure has more/less electron density separated by those vectors compared with the average structure. The types and signs of features can be directly related to the types of local correlations (16), and indeed the 3D-ΔPDF has been used to solve the local order in several disordered crystals (2,17,18).

Vacancy ordering
To remove the effect of atomic movements and isolate the effect of local vacancy ordering, the features in the 3D-ΔPDF can be integrated. As shown by Roth et al. (17), the integral amplitudes of features in the 3D-ΔPDF relate only to the substitutional disorder. The integrated peak amplitude of a peak at position ' is proportional to , where the summation runs over all pairs of atoms ( , ) separated by vector = ′, and is the difference in number of electrons of atom in the real structure compared to the average periodic structure.
It is possible to use this method on the current samples as the features in the 3D-ΔPDF are well-separated.
Details of the peak integration is shown in the supplementary materials.
The integral peak amplitudes are shown in Figure 3c and 3d. A positive value is found at the origin since atoms are separated by the zero vector to themselves. Surrounding the origin are strong negative integral amplitudes corresponding to nearest and next-nearest vectors in the Nb substructure (coordinates (½,½,0) and (1,0,0) ). This shows that there is a lower-than-average probability of finding two atoms separated by these vectors. This means that the structure tends to avoid nearest and next-nearest vacancy pairs. Then at slightly longer distances there are several positive amplitudes showing the preferred distances between vacancies. Again, the slowly cooled samples have strong correlations to much longer distances than the quenched samples.
The integral amplitudes can then be directly compared to the calculated 3D-ΔPDF obtained from a simulation of the vacancy distribution using the vacancy repulsion model (10). The model gives a positive energy to each nearest and next-nearest vacancy pair with the total energy of the system given by Here, 1 is the number of nearest-neighbour vacancy pairs and 2 is the number of next-nearest-neighbour

Structural relaxation around vacancies
The measured 3D-ΔPDF ( Fig. 3a and 3b) show strong features from local structural relaxation around vacancies. At (x,y,z) = (½,0,0), the 3D-ΔPDF is negative toward the center and positive away from the center.
This is the vector between Nb/vacancies and Sb. This means that when a Nb is present, the Sb will be slightly further away, and when a vacancy is present, the Sb will move towards the vacancy position.
In the average structure each Sb is surrounded by 6 close Nb/vacancy positions, forming an octahedron.
Avoidance of nearest and next-nearest vacancy requires each such octahedron to have at most one vacancy out of 6 corners. For x=1/6 there will be one and only one close vacancy per Sb in the ground state. This means that each Sb simply moves towards the one vacancy neighbour it has by approximately 0.144 Å.
Consequently, the (Nb,Sb) substructure can be seen as being built from Nb5Sb square pyramidal units, where the Sb is displaced slightly towards the empty Nb site, as illustrated in Fig. 1c.

Discussion
The quenched samples show more diffuse-like scattering than the slowly cooled samples, which have sharp additional peaks. This shows that the degree of local order can be tuned through the synthesis conditions (10). Surprisingly, the refined chemical compositions are virtually identical in all samples, x≈1/6, even though the nominal stoichiometry used in the syntheses differ substantially. Accurate determination of the composition was made possible by the high-quality single-crystal x-ray scattering data, avoiding errors from other impurity phases known to occur in these systems (9). This suggests that the stable composition is Nb5/6CoSb, and that there is very little room for variation in composition. Previous studies have argued that the ideal composition is Nb4/5CoSb based on simple valence electron counting (Zintl) rules (7). However, x = 1/6 is the highest vacancy concentration for which both nearest and next-nearest vacancy pairs can be avoided completely (10). This means that the energetic penalty for disturbing the local chemical bonding in the locally ordered vacancy structure is higher than the presumed gain in electronic energy expected from simplistic electron counting. It is noteworthy that the average structure is identical for all five samples, and also the local structural relaxation is the same with Sb is being off-centered by 0.144 Å and compositions deviate from 1/6. The measured transport data strongly reflect the presence of vacancies in their temperature behaviour, but presumably it is the lattice thermal conductivity, κL, which most directly is affected by the local structure. At room temperature, the samples differ by about 15% in κL and this could indeed be a direct effect of different vacancy distributions. It could therefore be advantageous attempt to use these differences in properties. Electron-diffraction data showed the diffuse features to be stable even at around 1000K, suggesting the vacancy distributions to be stable to quite high temperatures (9), making it possible to make devices with a specific vacancy order.
In summary, we have shown that defective half-Heusler materials Nb1-xCoSb have a strong tendency to vacancy ordering following a simple repulsion model. The optimal ordering of vacancies essentially fixes the stoichiometry of the samples irrespective of the nominal starting composition or synthesis method. Using analysis of both Bragg diffraction and diffuse scattering data with the 3D-∆PDF method, it is possible to quantify the structural relaxation around a vacancy site. Different local structure states are reached depending on the thermal treatment of the sample, and they appear to have appreciable effect on the thermoelectric properties. Advanced x-ray scattering techniques can unravel hidden local structures and for Nb1-xCoSb these local structures can be controlled by the synthesis conditions. If the local structure of crystalline materials can be more generally related to the properties, then a new frontier in materials research will be available. the other group of samples were slowly cooled using an induction furnace.

Tuneable local order in thermoelectric crystals
The quench samples are of the same type used in the published electron diffraction study (9). shows that only between 0.015 and 0.02 g are lost during synthesis of a 10 g sample, which corresponds to only between ⅓ and ½ of the added excess Sb.
Data was measured at the BL02B1 beamline at the SPring-8 synchrotron using a photon energy of 50.00 keV on a Huber four-circle (quarter chi) goniometer equipped with a Pilatus3 X 1M CdTe (P3) detector. For the quench samples a detector distance of 130 mm was used, while it was 260 mm for the slowly cooled samples to properly separate the sharp additional peaks.
Two measurements were made for each sample. First, a dataset for the strongest reflections were collected.
For this, a 600μm Ni film was used to attenuate the beam to 31% to avoid problems of too high flux on the detector. Here 4 runs were measured, each a 180° ω rotation with 900 frames, for χ=0° and χ =45° with the detector at 2θ=0° and 2θ =-25°. An exposure time of 0.8 seconds per frame was used. Then the weak scattering was measured without any beam attenuator. Here 6 runs were measured, each a 180° ω rotation with 900 frames, χ=0° and χ =45° with the detector at 2θ=0°, 2θ=-12.5° and 2θ =-25°. An exposure time of 1.6 seconds per frame was used. Finally the background and air-scattering was measured using the same set of exposure times, detector positions and beam attenuation as for the crystal measurement. For each combination of these, 200 frames of air scattering were measured.
To obtain the average structure, the images were converted to the Bruker .sfrm format (21) and the Bragg peaks were integrated using SAINT (22). The integrated data were processed and corrected using SADABS (23) and merged using SORTAV (24). Initial structure solution and refinement was carried out with SHELXS and SHELXL, using the Olex2 GUI (25,26), with subsequent structure refinement using Jana2006 (27) . The space-group is 4 � 3 . Anomalous scattering factors for 50keV were used, as implemented in Jana2006.
For the diffuse scattering analysis the data was converted to reciprocal space using a custom Matlab script.
During this process the data was corrected for polarization, the background scattering from air was subtracted, and a solid angle correction was applied. The resulting data was symmetrized using the 3 � For the production of a 3D-ΔPDF, the Bragg peaks of the average structure were punched and filled. Because the Bragg peaks and diffuse scattering do not overlap, the Bragg peaks were removed using a spherical punch based on the positions of allowed reflections of 4 � 3 . The punched holes were then filled using a 3D spline interpolation. The resulting data containing only the diffuse scattering/additional peaks was then Fourier transformed to give the 3D-ΔPDF.

Simulated data
Monte-Carlo simulations of the structure were performed using the Metropolis algorithm (19) using a custom python script. The model scattering was calculated using the Scatty software (28), and model 3D-ΔPDF obtained by Fourier transforming the model scattering using a custom python script.

Measured X-ray scattering
The measured X-ray scattering in the H0L and HHL planes are shown below for all samples. Figure S2: Measured X-ray scattering in the H0L and HHL plans for the different samples.
For the levitation-melt synthesis with nominal stoichiometry Nb0.84CoSb, data were measured on two crystals ("Q-0.84 # 1" and "Q-0.84 \# 2"). While crystal "Q-0.84 \# 2" shows quite broad diffuse scattering, crystal "Q-0.84 # 1" has sharper additional peaks. This suggests that the levitation-melt samples have inhomogeneities in the degree of short-range order, which might occur as different parts of the sample are cooled at different rates during the quenching process. In the vacancy repulsion model, vacancies on the Nb substructure avoid nearest and next-nearest pairs. For x = 1/6 there are several types of vacancy configurations that completely avoid nearest and next-nearest pairs (10). Of these, three have scattering patterns consistent with the measurement, as was shown previously (10). These are the A2, AD, and BD type ground states of the model. AD and BD are non-periodic and disordered.

Integrated peak amplitudes
The left panel of Fig. S5 shows the integrated amplitudes from the measured data on the SC-0.81 sample, which had the sharpest additional peaks in the scattering pattern. The three panels to the right show the calculated 3D-ΔPDF maps in the 010 plane for the three candidate ground states of the vacancy repulsion model without Sb and Co relaxations.
Both the A2 and BD models are in general good agreement with the experiment, while the AD model has several discrepancies in the sign of amplitudes. The black circles marks a feature which is different for the three models. For the BD model it is zero, for the A2 model it is positive and for the AD model it is negative.
In the experiment this integrated amplitude is approximately zero, but weakly positive, suggesting the "BD" model to be the most accurate. However, there are still some differences between the measurement and BD model, suggesting either the BD model is not the complete description, or that the peak integration method is not accurate enough. colour range compared to the third row to visualize the additional peaks better. The weakness of the peaks together with the following discussion suggest that the additional peaks in the experiment come from a different origin. The SC samples refine to stoichiometries of Nb0.822CoSb and Nb0.820CoSb, suggesting the deviation from = 1/6 to be quite small, if significant at all. With such small deviations the model does not produce clear peaks at the ring center. Furthermore, the peaks inside the 8-point rings in the measured data seem to become stronger with increasing q, while they decreases with q in the simulated data (see the blue, yellow and white circles). This could suggest that there is another reason for the peaks inside the rings. There are also further weak peaks in the measured scattering from the SC sample, not currently explained by the model, showing that there is still more to learn about these samples. Figure S6. Experimental data and calculated scattering for several simulated ground-state models.

Comments on the paper by Nan et al. (20)
A study by Nan et al. looks at high-angle annular dark-field imaging (HAADF), a type of high-resolution scanning transmission electron microscopy (STEM) measured on a sample of nominal composition Nb0.8CoSb (20). The HAADF is measured along the [110] zone axis where there is no blocking of elements of different kinds. Based on this it identifies that three different amplitudes correspond to the three different elements. It then tries to show that the short-range order is mainly displacive by making Fourier transforms of the measured and ideal amplitudes together with the measured and ideal positions of atoms.
The problem is that HAADF cannot see the vacancies. The transmission-type measurement gives a projection down along the 110 zone axis. This means that when a vacancy is present on a Nb position, the measurement will still see the Nb from other layers, as not all Nb down along a column will be missing. This is clear when looking at figure 2b in their paper. Clearly a dot is seen for every Nb position, as it is a projection along 110, giving an average over many layers. Because they cannot see the vacancies but only some degree of the displacements, they reach the conclusion: "We found that the short-range order is predominantly displacive." However, this is meaningless as their measurement does not allow analysis of the compositional disorder.
If the short-range order in the ED pattern was predominantly displacive, it would also not be compatible with the measured ED pattern ( fig. 1a in their paper). Displacive short-range order gives diffuse scattering which is very weak close to the beam center and increases away from the beam center, whereas compositional short-range order gives diffuse scattering which is more constant in intensity (for x-ray scattering it will decrease slowly away from the beamcenter because of the atomic scattering factor). This is illustrated in figure S7. The left column of fig. S7 shows the calculated diffuse scattering for the Nb/vacancy order, without any movements of Sb and Co. This gives diffuse scattering which decreases slowly in intensity due to the x-ray scattering factor. The central column shows diffuse scattering with only the Sb and Co displacements without any contribution from Nb. This gives diffuse scattering which is very weak close to the center and increases away from the center. The right column shows the full model with both the Nb/vacancy order and Sb and Co displacements, giving the best agreement with experiment.
In addition, the paper (20) avoids mentioning that the electron-diffraction data had already been explained to a higher accuracy using just substitutional short-range order (9), even though the authors were aware of this, as evident from them citing the paper in another context, and re-using several references that were first linked to the Nb1-xCoSb system in that paper.