research papers
Lower uncertainty bounds of diffractionbased nanoparticle sizes
^{a}Department of Applied Physics and Applied Mathematics, SEAS, Columbia University, 500 W 120th Street, New York, NY 10027, USA, and ^{b}Department of Mechanical Engineering, Özyeğin University, Çekmeköy, İstanbul 34794, Turkey
^{*}Correspondence email: icn2@columbia.edu
A selfconsistent analysis is reported of traditional diffractionbased particle size determination techniques applied to synthetic diffraction profiles generated with the Patterson approach. The results show that dimensions obtained from traditional techniques utilizing A_{L} versus L) plots, even when the entire rocking curve is used in the decomposition. It is recommended that, in routine analysis, all domain size determination techniques should be applied to all reflections in a diffraction pattern. If there is significant divergence among these results, the `average particle size(s)' obtained might not be reliable.
or Fourier analysis for singlecrystal nanoparticles have bestcase error bounds of around 5%. For arbitrarily shaped particles, lower error magnitudes are possible only if the zeroes of the thickness fringes are used. The errors for sizes obtained by integralbreadth and Fourierdecompositionbased techniques depend on the shape of the diffracting domains. In the case of integralbreadth analysis, crystal shapes which scatter more intensity into the central peak of the rocking curve have lower size errors. For Fourierdecomposition analysis, crystals which have nonuniform distributions of chord lengths exhibit nonlinearities in the initial ranges of the normalized Fourier cosine coefficient versus column length (Keywords: diffraction; particle size determination; Scherrer approach; integral breadth; Fourier analysis.
1. Introduction
Diffractionpeakshape analysis techniques for determining crystallite size (also termed domain size) are very popular, possess a massive literature base and, in some cases, have NIST Standard Reference Materials (SRMs) (Cline et al., 2020). However, the reliability, accuracy, precision and usefulness of the dimensions obtained from these techniques are still the subject of vigorous scientific debate. Some of the issues raised are common to all particle size determination methods. First, commonly used terms such as particle size, average particle size, particle size distribution etc. depend on the technique used for their determination and are not well defined for arbitrary particle shapes or characterization techniques (Matyi et al., 1987). Second, different experimental techniques used for size characterization are sensitive to different parts of a given particle size distribution and can yield vastly different values for the `average' size of the same assembly of particles. Thus, while none of these values are `wrong', none, individually, can adequately describe the size of a particle of arbitrary shape, or an aggregate of such particles of different sizes. Further, such differences usually become more pronounced if, in addition to a size distribution, there is a distribution of particle shapes within the aggregate. In such cases even computed distributions of size within the aggregates might be technique dependent.
Disagreements are also encountered when one compares crystallite sizes obtained from different diffraction formulations. For unstrained single crystals scattering in the kinematic regime, the crystallite size obtained from formulations based on the Scherrer approach is the maximum realspace length^{1} of the coherently scattering domains along the scattering wavevector, (Patterson, 1939a,b; Williamson & Hall, 1953). In contrast, Fourier analysis techniques of line shapes from such samples yield the average chord length along k (Guinier, 1994). For a given line profile, whether this average chord length obtained from Fourier methods is equal to the maximum crystal dimension obtained from the Scherrer approach depends on the shape of the crystal.
If there are additional sources of broadening in the diffraction line profile, as well as contributions from neighbouring peaks and the background intensity profile, obtaining an accurate and useful crystallite size from diffraction line profile analysis becomes more laborious and more error prone. This may be one of the reasons for the many conflicting statements in the literature about the efficacy of various formulations. For example, while some authors warn against the use (abuse) of the Scherrer equation (Matyi et al., 1987; Palosz et al., 2010; Leoni, 2019), others have concluded that it is a reliable technique, especially if strain broadening is negligible and lowerorder peaks are used (Kaszkur, 2006; Ying et al., 2009; Dorofeev et al., 2012). Some authors discount the use of diffraction techniques altogether for specific applications. Tomaszewski (2013) asserted that critical grain sizes for sizeinduced phase transitions cannot be determined using diffraction methods, stating that `none of the known methods give the correct value for the crystallite size of nanocrystals! It is only possible to talk about the range of values.'
Another factor which complicates the comparison of crystal sizes obtained from various diffraction formulations is the lack of rigorous uncertainty values. Many experimental parameters influence the peak profile, some in a highly nonlinear manner, and the propagation of errors for some formulations is quite complicated. Currently there are no rigorous formalisms which can compute, from first principles, the full uncertainty range associated with crystallite sizes determined by diffraction methods (Young et al., 1967).
In this work we have investigated the lower error bound (minimum uncertainty) associated with the accuracy and trueness (closeness between the average of an infinite number of replicate measured quantity values and a reference quantity value; Noyan et al., 2020) of average crystal sizes obtained from the most popular diffraction formulations. For these computations we utilized the Patterson equations (Patterson, 1939a,b) to generate standalone diffraction peaks for various reflections of perfect Au crystallites with thinfilm, cylindrical and spherical geometries. These intensity data are exact, with no counting uncertainty; thus, the peaks are exactly reproducible with infinite precision. These peak profiles were analysed for the average `particle size' using the Scherrer equation (both FWHM and integralbreadth modes) and Fourierdecomposition approaches. Also, since the Patterson equations are based on exact geometric shapes, and the computed peaks are free from all confounding effects (such as limited angular range, background intensity, peak overlap and asymmetry), we expected this approach to provide the best achievable accuracy, and hence the minimum error bound, for particle sizes obtained by diffraction techniques.
2. Theory and simulations
We first reprise the 1939 Patterson formulation of the Scherrer equation (Patterson, 1939a,b) and compare its predictions of particle size with those of the Fourier analysis of line shapes (even though most of this material is available in the literature, it is scattered across many publications which are, mostly, somewhat terse treatments).
The Patterson formulation considers only kinematic diffraction which neglects multiwave scattering and assumes that the energy diverted into the diffracted beam is negligible. Under these assumptions, for a plane wave incident on a crystal, the wave scattered by an atom at position will have a phase difference from the wave (with the same wavevector) scattered by an atom at the origin. The total diffracted wave amplitude in the direction of is the sum of the waves scattered by all the atoms in the crystal. Since Xrays are scattered by the core electrons surrounding the atoms, it is more accurate to consider the electrondensity distribution in the crystal rather than the atomic positions,
Here, is the envelope function of the sample crystal, which is equal to 1 inside the crystal boundary and zero outside, and is the electrondensity distribution of a triply periodic infinite crystal,
v and F_{hkl} are the volume and respectively, of the of the crystal, and is the reciprocalspace vector.
The total diffracted wave amplitude is the sum of diffracted waves from all electrons over the entire crystal, taking into account the proper phase differences. This can be written as
Here, is the wavevector, and ), termed the shape function, , is the Fourier transform of the sample envelope function ,
Equation (3) indicates that the shape of the Bragg peak in angular space is determined only by the realspace shape of the diffracting crystal. The type and chemistry of the modulate the amplitude at a given reciprocalspace position (or diffraction angle) mainly through the F_{hkl}.
If the crystal is infinitely large in all directions, equation (3) yields Bragg's law; in this case and finite amplitude (and thus intensity) will only be observed at the exact Bragg condition where .
For finitesized crystals the shape of the diffraction peak, i.e. the angular distribution of scattered intensity around the must be obtained by evaluating the triple integral in equation (4). Patterson and many others have computed for many regular geometric shapes. In this study we will only consider three simple cases where the diffracting crystal is (i) a thin film, (ii) a cylinder and (iii) a sphere, respectively. The corresponding geometries are sketched in Fig. 1.
2.1. Diffraction from a singlecrystal thinfilm slab
We consider a radial scan from a perfect singlecrystal slab of (realspace) thickness t_{f} and infinite inplane area, where the diffracting (hkl) planes are parallel to the surface of the film with their normal [hkl] coincident with the normal to the film surface [Figs. 1(a) and 1(b)]. In this case evaluation of equation (4) yields the normalized cardinal sine (sinc) function,
where V is the volume within the singlecrystal thin film illuminated by the Xray beam. By substituting this expression into equation (3) we obtain the amplitude at a given . The corresponding diffracted intensity along is obtained by multiplying the amplitude by its conjugate. Expressing in terms of deviation from the , we obtain
Here, N is the number of unit cells along and k = 1/λ. Equation (6) shows that the shape of the (kinematic) diffraction peak from a perfect singlecrystal thin film illuminated by a plane wave is described by the squared normalized sinc function^{2} = . When the thickness t_{f} approaches infinity, the intensity will be zero at all Δ2θ except Δ2θ = 0. For a finite thickness, the function will have a primary peak with its maximum intensity (I_{max} = 1) at Δ2θ = 0, which is bracketed by symmetric subsidiary peaks (thickness fringes) of diminishing intensity with increasing Δ2θ. An example pattern computed for a 5 nm thick Au film using Cr Kα radiation is shown in Fig. 2. One can use three techniques to obtain the film thickness from this peak profile using equation (6):
(i) The periodicity of the thickness fringe zeroes, T, is related to the film thickness^{3} by
(ii) The halfmaximum intensities of the diffraction peak will be reached when = ±1.39. Thus, the FWHM β of the primary peak will given by
This is the Scherrer equation for FWHM, with the Scherrer constant C = 0.885.
(iii) The integral breadth β_{I} of the function is related to the film thickness by
This is the integralbreadth form of the Scherrer equation with a Scherrer constant of unity.
While the forms of equations (7a)–(7c) are very similar, their efficacies in computing the film thickness are not. As we showed previously using highresolution diffraction from silicononinsulator films (Ying et al., 2009), equation (7a) is the easiest to use, and the most accurate, for singlecrystal thin films since it does not utilize a `shapedependent' constant and the zeroes of the thickness fringes can be determined very precisely when highresolution XRD systems on highintensity sources are utilized. The FWHM form of the Scherrer equation requires the Scherrer constant for the particular geometry. Also, the FWHM obtained from an approximate function, such as a Gaussian, fitted to introduces a small error in β and hence in t_{f,β}.
Equation (7c) also does not require a Scherrer constant. However, computation of the integral breadth β_{I} is nontrivial. Fitting only a portion of the profile with an approximate function, such as a Gaussian (dashed line in Fig. 2), will cause errors in the computed film thickness t_{f,IB} since (i) the tails of the Gaussian function do not fit the tails of the central peak of the function and (ii) fitting only the central peak underestimates the integrated intensity of the radial scan; this would result in a film thickness t_{f,IB} larger than its actual value.
To investigate this effect we computed β_{I} by evaluating the integral
over various – ranges in Fig. 2. For ease of interpretation we chose these ranges to include the central peak plus identical numbers of satellite peaks, bounded by their respective zeroes, on each side. In Fig. 3 we plot the t_{f,IB} values corresponding to these different integration ranges, starting with the range for just the central peak and ending with the range with the central peak plus eight satellite peaks on each side. In this figure we also include the fractional excluded ΔP_{A} for each such range,
If just the integral breadth of the central peak is used to compute the film thickness, t_{f,IB} is ∼10% larger than the actual film thickness, since the area of the central peak is ∼10% smaller than the true total integrated intensity, . As the angular range of integration increases, I_{INT} and t_{f,IB} asymptotically approach their ideal values; analysis of the range containing the central peak plus four satellites at each side (sixth point from the left) yields ∼2% deviation.
To investigate this issue further we fitted multiple Gaussians to the profile shown in Fig. 2. An example with seven peaks (central peak plus three satellites on each side) is shown in Fig. 4. The t_{f,IB} values obtained using the sum of the integral intensities of all fitted peaks are summarized in Table 1. For ease of comparison the corresponding t_{f,IB} values from the integration of equivalent ranges of the profile are also listed. We observe reasonable agreement between the two approaches.

Table 1 also shows that the film thickness values obtained from the FWHM form of the Scherrer equation, equation (7b), do not change significantly with increasing number of fringe peaks fitted simultaneously with the central peak. All of these values are within 5% of the actual (ideal^{4}) thickness. As shown in the last row of Table 1, one obtains the best accuracy when using the zeroes of the thickness fringes [equation (7a)] to compute the film thickness t_{f}.
To check the generality of these observations, we computed the Patterson profiles for the 111 peaks of singlecrystal thin films of various thicknesses and repeated the computations described above for each profile. The results are summarized in Fig. 5, where the fractional thickness error Δt = (t_{f,IB} − t_{f})/t_{f} and the corresponding ΔP_{A} values [equation (9)] are plotted versus the normalized integration range for all cases. Both of these parameters are independent of size. In all cases the fractional error obtained from using the integral breadth obtained by fitting a single Gaussian to the central peak yielded an identical fractional error corresponding to the first Δt_{f,IB} datum (integration of the central peak). For all films, when t_{f,IB} was computed using the definite integral , the result had zero fractional error.
We conclude that it is challenging to obtain the true integral breadth from fitting a single function of any approximate form to the central peak of a radial scan. Numerical integration of an experimental profile over the largest possible angular range might yield better results. However, given the asymptotic approach of the integrated intensities to the definite integral values in Fig. 5, errors of several percent would be expected in the computed thickness values, since the presence of neighbouring peaks or experimental issues might limit the angular scan range.
2.2. Diffraction from singlecrystal cylinders and spheres
In the case of a cylinder of height L_{cy} and base radius R_{cy}, where is along the radial direction, equation (4) becomes
The corresponding diffracted intensity in terms of deviation from the Δ2θ is
The halfmaximum intensity is reached when = ±0.808. The FWHM (β) value can be expressed in terms of the base diameter D_{cy} of the cylinder as
Likewise, for diffraction from a singlecrystal sphere of radius R_{sp} (diameter D_{sp}), we obtain
Equations (11) and (13) yield radial scan profiles similar to the thinfilm profiles [equation (6)]. In Fig. 6 we plot the 220 radial scans computed with λ = 2.2909 Å for a cylinder and a sphere, with D_{cy} = D_{sp} = 5 nm, and for a thin film, with t_{f} = 5 nm. All three profiles have a central (primary) peak of unit amplitude bracketed with symmetric thickness fringes. The FWHMs of the central peaks are slightly different for the three sample geometries, commensurate with the differences in the Scherrer coefficients in equations (7b), (12) and (14). The rates of decay of the thickness fringes are significantly different, with those belonging to the sphere profile extinguishing fastest. The periods of the zeroes of the thickness fringes are equal for all three cases, i.e. the thicknessfringe zeroes are independent of shape. Consequently, equation (7a) can be used directly for computing D_{cy} and D_{sp}.
In contrast to the thinfilm case [equation (6)], the integral breadths of the peaks described by equations (11) and (13) cannot be used directly in equation (7c) to calculate D_{cy} and D_{sp} since these equations are not normalized (while the intensities generated by both functions are in the range between 0 and 1, their norms are not unity). In the case of a cylinder, the integrated intensity of the rocking curve is given by
Thus, the integral breadth obtained from the rocking curve of the cylindrical sample is related to the base diameter of the cylinder by
For the case of a spherical sample, a similar treatment yields
Consequently, when the Patterson equations are used to simulate radial scan profiles from spherical or cylindrical singlecrystal samples, a `Scherrerlike' proportionality constant (shape constant) must be included when linking the integral breadth of these profiles to the diameters of such samples, D_{IB,cy} and D_{IB,sp}, even when the radial scans are integrated over infinite angular ranges, .
Integration over a narrower angular range^{5} will introduce a truncation error in the diameters D_{IB} even when such a proportionality constant is used (Table 2). For a given normalized integration range this error will be smaller for cylindrical and spherical single crystals than singlecrystal thin films, since the central peak of these patterns contains a larger fraction of the diffracted intensity (Fig. 7). For both of these geometries, evaluating the diameter from a Gaussian fit with the theoretical Scherrer constants to the central peak only yielded ∼2% error (Table 2). This small error is due to the Gaussian peak approximation. If we calculate the Scherrer coefficients by correlating the FWHMs of the simulated peaks fitted with Gaussians to the actual sphere diameters used in the simulations, we obtain C = 1.12 and 0.99 for the sphere and cylinder coefficients, respectively. The corresponding theoretical values for these geometries obtained from the exact functional forms are C = 1.16 and 1.03 (Patterson, 1939a,b).

3. Fourier analysis
Fourier decomposition is also widely used to obtain crystallite size, lattice strain and ; Warren, 1969; Klug & Alexander, 1974; Schwartz & Cohen, 1987; Guinier, 1994) and in journal articles (Warren & Averbach, 1950, 1952; Langford et al., 2000; Ida et al., 2003; Lucks et al., 2004; Scardi & Leoni, 2006; Mittemeijer & Welzel, 2008; Cline et al., 2020). The basic approach consists of correcting the measured diffraction profiles of one or more diffraction peaks for background, peak overlap and instrumental effects. These corrected profiles are then decomposed into their Fourier components, from which the relevant parameters are computed. In the absence of strain (or other samplerelated) broadening, the Fourier cosine coefficients of any reflection can be used to determine either the distribution of crystallite sizes or an average crystallite size for a polycrystalline powder. For a powder sample with unimodal size and shape distribution, or for a kinematically scattering single crystal, the average crystallite size determined from Fourier decomposition is the average chord size along the diffraction vector , which depends on the realspace shape of the crystallite, the diffraction geometry and the of the reflection used in the analysis (Guinier, 1994). In what follows we apply the basic Fourier analysis directly to synthetic peak profiles generated using the Patterson equations for thinfilm, spherical and cylindrical samples. Rather than reprising the theory, which is readily accessible in the literature, we will provide a comprehensive discussion of the actual steps of the basic analysis which yields the average dimension along , with further details about the use of the Fourier transform provided in Appendix A.
distributions from diffraction peak profiles. The fundamental theory and its extensions are widely discussed in textbooks (Taylor, 1961The discrete peak profiles used in the analysis were computed using a simple Mathematica notebook (Wolfram, 2021), which output the simulated diffraction peak intensities as [I_{j}, Δ2θ_{j}] data sets with M equidistant δ2θ points. Here, I_{j} is the jth intensity point at angular position Δ2θ_{j} = j*δ2θ, where . We then used a discrete Fourier transform (DFT) formulation restricted to real numbers, , to express each peak as a periodic function I_{F}(Δ2θ) in the angular range (0, Mδ2θ),
Here, and are, respectively, the Fourier cosine and sine coefficients of order , and I_{F}(Δ2θ_{j}) is the (Fourier) synthesized intensity at angular position Δ2θ_{j}. While the profiles generated by the Patterson equations are centred at Δ2θ = 0, the peak profiles obtained from equation (19) are shifted by half of the angular range and have peak maxima at Δ2θ = (Mδ2θ)/2. This shift does not affect the analysis results.
In Fig. 8(a), the cosine and sine coefficients and computed for the 111 reflection of a hypothetical 5 nm thick Au thin film using Cr Kα radiation are shown for n ≤ 452. The 111 peak profile was sampled at 905 steps (M = 905), 131 of which were in the central peak. For all n, , reflecting the symmetry of the Bragg peak. The function I_{F}(Δ2θ) computed using these coefficients showed good fidelity to the actual [I, δ2θ] data set [Fig. 8(b)]. The maximum residual normalized intensity difference between the actual and recalculated intensity data sets was 5 × 10^{−7}.
To obtain the thinfilm thickness, we computed the column length L (the length of the columns of the unit cells in the direction of the diffraction vector; Guinier, 1994) as for each Fourier order n, and the corresponding absolute values of the normalized cosine coefficients,^{6} . Here, is the full range for [equation (6)] over which the Patterson intensities were computed. In what follows we will omit the order subscript on the column length for brevity, following the common usage in the literature. Fig. 9 depicts the variation in A_{L} versus L for the full set of thinfilm simulations. We observe that
(i) for all film thicknesses, A_{L} versus L varies linearly for L < t_{f}, and
(ii) leastsquares lines fitted to A_{L} versus L data over the entire 0 < L < t_{f} range, or any of its subsets, intersect the abscissa at L = t_{f}.
These observations agree with our experimental results obtained from singlecrystal silicononinsulator thin films (Ying et al., 2009).
In contrast to the thinfilm plots shown in Fig. 9, A_{L} versus L plots for spheres and cylinders exhibit sigmoidlike profiles (Fig. 10) with quasilinear central regions. For a cylinder simulation with D_{cy} = 40 nm, the linear A_{L} versus L range is approximately 46% of the linear range for the 40 nm film. For a sphere simulation with D_{sp} = 40 nm, the linear range drops to 30%. Thus, the geometry of the scattering volume (or, more precisely, the form of the distribution of chord lengths along the diffraction vector within the crystallite volume) introduces nonlinearities into the initial and final ranges of the A_{L} versus L plots. Consequently, the curvature in the initial part of such plots cannot be ascribed solely to the traditional `hook effect' (Warren, 1969) caused by background subtraction issues or incomplete measurement ranges.
The presence of these nonlinear regions requires the use of a consistent technique for defining the linear leastsquares fit range to ensure precise determination of the desired crystal dimension. For this purpose we computed the (local) slope of the line segments for each adjacent pair of normalized cosine coefficients (A_{L}_{j}, A_{L}_{j−1}), and in the regression analysis we included only those pairs yielding local slopes within 10% of the mean slope of the line fitted to the middle `linear' region. These regions are highlighted in Fig. 10. The lines fitted to these regions intersected the abscissa at L values quite close to the average chord lengths of the cylinder and sphere shapes along the diffraction vector (Guinier, 1994),
For symmetric radial scans obtained from singlecrystal thin films, the average cord length is identically equal to the film thickness,
To obtain better statistics we applied the Fourier formalism to spheres and cylinders with diameters ranging from 5 to 40 nm and obtained the Fourieraveraged cylinder and sphere diameters, and , respectively, from the intercepts of the linear regression lines. As shown in Figs. 11(a) and 11(b), both and depend linearly on the respective (geometric) diameters, with and . In both cases there is a ∼4% deviation from the theoretical values predicted by equations (20a) and (20b) for the ranges analysed. This error is larger than the `fit' errors assigned to the slopes by regression analysis.
We were not able to obtain satisfactory diameter values when we used the initial regions of the A_{L} versus L plots for these geometries as suggested in the literature. We note that selection or exclusion of a few points at the start and/or end of the linear range could change the slope by an additional 5% or more, with concomitant changes in the and values. Consequently, for these specimen geometries, even if all other broadening sources were successfully eliminated and there was no surface roughness, experimentally determined particle sizes could have uncertainties in the region of ±5%.
3.1. Multiple peak comparison
As a final test of the Fourier analysis formalism we investigated the effect of the curved A_{L} versus L regions in two orders of the same reflection, postulating that, if our approach to determining the `linear region' is correct, both reflections should yield the same average chord length within the `fit' error since there is no strain broadening. Fig. 12(a) shows the 111 and 222 reflections of a 40 nm Au singlecrystal sphere computed using equation (13) with λ = 2.2909 Å. Due to the cosθ_{B} term, the FWHMs of the primary peaks and the period T of the fringe minima are different. These differences are reflected in the (raw) cosine and sine coefficients, and , of the Fourier series corresponding to these peaks [Fig. 12(b)].
In contrast to Fig. 12(b), the normalized cosine parameters A_{L} for the two reflections are very close in value and exhibit similar L dependencies (Fig. 13). Also, the local slopes of the two curves are almost identical. Consequently, the effective sphere diameters [equal to the mean chord length of the sphere, equation (20b)] obtained from extrapolation of the linear portions of these A_{n} versus L plots are very close for the two reflections (Table 3), reflecting the absence of (orderdependent) strain broadening in our model.

4. Summary
In this study we first simulated standalone diffraction peaks from kinematically scattering strainfree single crystals in the shape of thin films, cylinders and spheres using their Patterson equations. These simulated peak profiles were centred at the exact β) and integralbreadth forms of the Scherrer equation, (ii) the period of the thickness fringes, and (iii) singlepeak Fourier analysis. The difference between the dimension input into the simulation and that obtained from a given technique was taken as a measure of the minimum uncertainty (error) associated with the results from the particular technique. To obtain better statistics this analysis was conducted for a range of crystal sizes. We observed the following:
with tails extending to , and were free from other samplebased broadening effects such as dislocations, stacking faults, strain, instrumental broadening or background profiles. We then computed the relevant crystal dimension, thickness or diameter, along the scattering vector for these crystals using (i) the FWHM ((i) The best accuracy was achieved when the zeroes of the thickness fringes were used to compute the crystal dimension along . For well defined fringes, subångström precision in the relevant dimension could be achieved. This observation agrees with our previous experimental work on silicononinsulator singlecrystal thin films.
(ii) The use of a Gaussian function to approximate the central peak of the diffraction profile introduced errors into the crystal dimensions obtained from both forms, FWHM and IB, of the Scherrer equation.
(a) Use of the Gaussian FWHM in the Scherrer equation yielded approximately 5% error in the computed crystal size. This error was independent of shape.
(b) A shape constant, not equal to one, was needed in the computation of the diameters of cylindrical and spherical crystals using integralbreadth values, since the Patterson equations describing the peak profiles for these shapes are not normalized functions, even though the computed intensities range between zero and unity.
(c) Even when such a shape constant was used, if the relevant crystal size was computed from the integral breadth of a Gaussian function fitted to the central peak, errors in the range of 10% (for thin films) to 2% (for spheres) were observed. These errors are due to the approximation of the diffraction line shape with the Gaussian function and depend on the shape of the crystal. Similar errors are possible when other approximation functions (such as pseudoVoigt, Pearson VII or Lorentzian formulations) are employed.
(iii) Singlepeak Fourier techniques also yielded shapedependent crystal size errors.
(a) For perfect thin films, the variation in A_{L} versus L was linear for all thinfilm dimensions used in the study (Fig. 9). Average thickness values computed from the intercepts of these plots on the abscissa had negligible errors.
(b) A_{L} versus L plots computed for cylindrical and spherical specimens had sigmoidlike A_{L} versus L profiles (Fig. 10).
(1) The initial concavedown curvature in these plots was caused by the shape of the crystal. It is not due to the traditional `hook' effect and could not be corrected using traditional approaches. The Fourier transform of a diffraction peak from a crystallite yields the autocorrelation function of the crystallite along the scattering vector. This is analogous to the distribution of signal energy across frequencies in signal processing and is proportional to the average number of unit cells across the crystallite Δ2θ). This number is constant when the crystal is a thin film. For spherical or cylindrical samples, the variation in this number with Δ2θ is never truly linear. Thus, the shape of the A_{L} versus L plots will depend on the shape of the crystal; if the shape of the crystal is nonsymmetric along the scattering vector with Δ2θ, the A_{L} versus L plot will also be nonsymmetric.
participating in diffraction for the particular scattering angle (expressed as the deviation from the exact(2) The presence of these nonlinear regions necessitated fitting the central region of these plots with a line to obtain the relevant crystal diameters. Since the definition of these linear regions was somewhat arbitrary, selection or omission of a few points could contribute errors of about 5% in the extrapolated average crystal diameter.
(iv) As shown in Appendix A, when Fourierdecomposition programs which are optimized to treat real number data are applied to a peak profile with M intensity points, the cosine and sine coefficients for 1 ≤ n < M/2 will be twice the magnitude of the coefficients obtained from standard programs, while the coefficients for the zeroth order and the Nyquist frequency (n = 0, n = M/2) will be identical. For such codes the zerothorder (raw) cosine coefficient must be corrected if it is to be included in the analysis.
(v) For Fourier analysis, all other considerations being equal, it was much better to use a broad angular range with fewer steps, rather than a narrow angular range with finer steps. If a narrow range is used, using `padding' might enable the recovery of useful data from the analysis (Appendix A).
5. Conclusions
Even for the best case, where all confounding issues have been eliminated, classical diffractionbased particle size values, with the exception of those from thickness zeroes, have fractional errors of around 5% for the coherent domain size along the scattering direction. This value is a lower bound; the presence of other linebroadening sources, instrumental effects or emergent scattering artefacts such as the nanoparticle size error (Xiong et al., 2018, 2019; Kaszkur, 2019) will cause larger uncertainties. On the basis of our results we recommend that, once a diffraction pattern has been corrected for instrumental broadening, all domain size formalisms, i.e. thickness fringe zeroes,^{7} FWHM and IB formulations of the Scherrer equation, and Fourier decomposition, should be applied to all available reflections; these results should then be evaluated together. If all values agree within experimental error, then the sample probably exhibits a uniform chord length distribution along the scattering vectors. If there is significant divergence among these results, the `average particle size(s)' obtained might not be reliable. In such cases, if other sources of broadening can be eliminated, more sophisticated approaches involving domain size distributions such as those described by Scardi & Leoni (2006) might be undertaken.

APPENDIX A
Fourier decomposition of Bragg peaks
This analysis starts with the measurement (or computation) of a given peak profile spanning an angular range Δ2θ_{FR}. In most literature treating Fourier analysis of signals, this is termed a `frame'. We assume that the scattered intensities in this frame are sampled at M equidistant angular positions, yielding an array of intensity versus angular positions, [I_{j}, Δ2θ_{j}], where Δ2θ_{j} = jδ2θ, , and Δ2θ_{FR} = Mδ2θ. This intensity profile is corrected for instrumental effects and background as needed, and then the intensity values are normalized such that . The first goal of the analysis is the determination of a Fourier interpolation function I_{F}(Δ2θ), which matches the (normalized) intensities for all Δ2θ_{j},
where n is the order of the Fourier coefficients. Since at least two data points are required per period to resolve a wave, the highest unique `frequency' component resolvable in the Fourier transform of the diffraction pattern is at n_{NL} = M/2, which is termed the Nyquist limit (Smith, 1999; Avillez et al., 2018).^{8}
The cosine and sine Fourier coefficients are related to the complex coefficients through
The complex Fourier coefficients can be computed from
The corresponding cosine and sine coefficients can then be computed from equation (22). In the literature, equations (21) and (23) are termed the discrete Fourier transform (DFT) and the inverse discrete Fourier transform, respectively.
A1. Decomposition algorithms
In practice, freely available or commercial programs are used to compute cosine and sine coefficients. Most of these programs have several options for data conditioning and use different algorithms. Consequently, the number and magnitude of coefficients obtained from different programs may not agree. We compared three programs for this purpose.
In the first case, since for all j, we utilized a DFT code (Project Nayuki, https://www.nayuki.io/page/freesmallfftinmultiplelanguages) limited to data in the realnumber domain. For such data both cosine and sine coefficients are real and . Consequently, = and = for 1 ≤ n < M/2, and all coefficients are scaled by the number of sampled intensities M. Using these coefficients one can compute the original Bragg peak from the inverse DFT equation:
Here the summation is up to the Nyquist limit, n_{NL} = M/2.
In Fig. 14 the 111 peak profile obtained from the Patterson equation for a hypothetical 50 Å thick Au thin film using Cr Kα radiation is shown, with Δ2θ_{FR} = 0.725 rad (approximately ∓8β) and δ2θ = 0.0008 rad; this yields M = 905 intensity values for the DFT analysis.
Computing the Fourier coefficients up to the Nyquist frequency n_{NL}, we obtain 453 cosine and sine coefficients, which are plotted in Fig. 15(a). The sine coefficients are negligible for all n, reflecting the symmetry of the peak profile. The cosine coefficients vary linearly with n for 1 ≤ n ≤ 13 and decay to negligible values for n ≥ 14 [Fig. 15(b)]. The magnitude of the zerothorder coefficient , reflecting the integrated intensity of the peak profile over the measurement interval, is smaller than the absolute magnitudes of for 7 ≤ n. As we mentioned above, this is a direct result of a DFT program limited to data in the realnumber domain; when such programs are utilized in the analysis of experimental diffraction data for size and rootmeansquare strain analysis, one can either exclude from the set of Fourier cosine coefficients or multiply this value by two. In our work we exclude . If we substitute all cosine and sine coefficients shown in Fig. 15(a) in the inverse DFT equation [equation (24)] to resynthesize the 111 peak, we obtain an excellent fit with the peak profile computed using the Patterson equation (Fig. 14, red dots). In this frame the maximum intensity difference between the two traces is 4 × 10^{−7}.
In our second approach we coded the following fast Fourier transform (FFT) equations in MATLAB (MathWorks, 2019):
(i) For Fourier decomposition,
(ii) For Fourier synthesis,
Equations (25a) and (25b) yield M cosine and sine coefficients, M/2 of which are independent. Thus, when applied to the Bragg peak shown in Fig. 14, we obtain 905 coefficients each for and [Fig. 16(a)]. As expected, with the exception of n = 0 and n = M/2, the magnitudes of the cosine coefficients obtained from equation (25a) are one half of those obtained from equation (23) [Fig. (16b)]. This figure also intimates that obtained from equation (25a) will exhibit the traditional monotonic decline (linear decline in the case of thinfilm single crystals) from n = 0.
For this approach, all M coefficients computed from equations (25a) and (25b) are required to resynthesize the diffraction profile. Such generalpurpose programs consume more resources than programs optimized for real numbers.
Finally, we tested the FFT algorithms in the MATLAB and OriginPro (OriginLab, 2021) software packages. Both programs use algorithms based on the FFTW subroutines (http://www.fftw.org) and yield `unscaled' Fourier coefficients. These must be scaled (normalized) by the number of data points (i.e. multiplied by 1/M) before they can be used in resynthesizing the Bragg peak profile. (While scaled coefficients must be used for synthesizing peak profiles, coefficients from all programs can be used directly as input for Fourierbased average size analysis, since must be normalized (divided) by or for this purpose.)
A2. Specification of frame range and number of data points
Equations (21a)–(26) show that, since Δ2θ_{FR} = Mδ2θ, only two parameters of the DFT/FFT analysis can be independently specified to analyse Bragg peak profiles. In a recent article, de Avillez et al. (2018) investigated the optimization of these parameters for determining crystallite size distributions when confounding instrumental effects were present. Those authors did not extend the treatment to the simple Fourieraveraged sizes when instrumental effects were absent. Here, we use numerical analysis of the Bragg peak shown in Fig. 14 to investigate the influence of Δ2θ_{FR} and δ2θ on the accuracy of the average crystal dimension along the diffraction vector .
To determine the optimum Δ2θ_{FR}, δ2θ combination at minimum experiment and/or processing time, one can fix the frame range Δ2θ_{FR} and vary δ2θ, and then fix the optimum δ2θ obtained in the first step and change the frame range. We started by setting Δ2θ_{FR} = ∓8β, as shown in Fig. 14, and changing the step size. Relevant parameters are listed in Table 4. The average Fourier film thickness for each case (tabulated in the last column) was obtained from the intersection of the line fitted to its A_{L} versus L data. These are shown in Figs. 17(a) and 17(b). For all reasonable step sizes the average film thickness is within 5% of the true film thickness. Even when the step size was not reasonable – there were only five points in the central peak for M = 33 – the average Fourier thickness of the film was within 5% of its true realspace value. The only effect of increasing the step size is the introduction of some nonlinearity for large L, L < t_{f}. We conclude that 458 steps, corresponding to a step size 3.6% of the FWHM of the primary peak, are good enough for this ideal analysis.
For the second step of the optimization analysis we fixed the step size at 0.90° (0.00160 rad) and investigated the effect of frame range on . Relevant details for this analysis are summarized in Table 5. We observe that decreasing the frame range dramatically curtails the number of Fourier coefficients N_{L}^{*} which yield column lengths in the linear range for the A_{L} versus L plots (for n ≥ 1). Even when the entire central peak is sampled with 66 steps (second row of Table 5) we cannot perform the analysis. On the other hand, if the range includes just the central peak plus one pair of (bracketing) thickness fringes, the extrapolated value is very accurate, even though the number of Fourier coefficients is quite small [Figs. 18(a) and 18(b)]. These observations imply that capturing most of the Bragg peak is quite important: for the same experimental duration it is better to use a larger step size (with appropriate counting time) over a broader frame range than to perform a finestep scan over a narrower angular frame.

Finally, we investigated the effect of `zero padding' on the analysis. In this procedure the frame range is artificially broadened by adding zero intensities to equidistant angular positions outside the measurement range. In our analyses we chose to extend the frames of simulations described by rows 1–4 in Table 5 to the full range shown in Fig. 14. Consequently, the Fourier transformation (DFT) acted on 458 intensity points in all cases. However, only a fraction of these intensity points centred on the primary peak were finite and `real'. The rest were `padding'.
`Padded' profiles containing the top 75 and 100% of the primary peak, with 40 and 66 original intensity values, M_{o} = 40 and 66, respectively, are shown in Figs. 19(a) and 19(b). The corresponding Fourier coefficients and the A_{L} versus L plots (for n ≥ 1) obtained from these coefficients are shown in Fig. 20. For both cases, artificially expanding the frame range by tacking on zero intensities resulted in the recovery of 14 finite A_{L} values for L ≤ t_{f}. However, Figs. 20(a) and 20(b) show that padding introduces oscillations in the values. These oscillations decay with increasing `real' frame range and are due to the artefacts introduced into the peak profiles by the padding. We also observe both oscillations and curvature in the A_{L} versus L plots (for n ≥ 1). These also diminished with increasing true frame range, vanishing completely for adequate sampling: for M_{o} ≥ 264, the A_{L} versus L plots with 454 real intensity points, and with 264 or more real intensity points padded to 454 total points, were coincident.
In Table 6 we list the average film thicknesses obtained from the A_{L} versus L plots (for n ≥ 1) by fitting the linear central ranges. The computed values are sufficiently accurate for most structural characterization applications. We conclude that Fourier analysis, properly implemented, is quite robust for the determination of average domain size in singlecrystal (or monodisperse) polycrystalline samples with no strain broadening effects, even when padding is used. On the other hand, ascribing nonlinearities in Fourier coefficient plots and their extension to structural features of the sample appears to be nontrivial. Such effects can also be caused by the type of algorithm used for computing the raw Fourier coefficients, and , or be due to asymmetries in the peak shape caused by frame range selection, padding or other nonsamplerelated reasons.

Footnotes
^{1}Some publications (Matyi et al., 1987; Leoni, 2019) identify this dimension as a `reciprocalspace length'. This issue is discussed later in the article.
^{2}In what follows we neglect the multiplier terms in the Patterson intensity equations for brevity. Consequently, the of the particular crystal will be represented in the peak shape only by the θ_{B} of the particular reflection being simulated.
^{3}In equations (7a)–(7c) we differentiate between the thicknesses obtained by these techniques using subscripts. In the absence of errors all of these terms must be equal to the realspace film thickness: t_{f,z} = t_{f,β} = t_{f,IB} ≡ t_{f}.
^{4}The Patterson formulations are based on ideal geometric particle shapes with smooth surfaces. The presence of surface roughness will exacerbate such errors and increase the uncertainty associated with the trueness of such values.
^{5}Since equation (17) is not well behaved in the vicinity of Δ2θ = 0, numerical integration of the function in this range requires care.
^{6}In general, A_{n} (the normalized cosine coefficient of order n) is not equal to the (raw) cosine coefficient . In our calculations we used for normalization (instead of the zerothorder coefficient normally employed for this purpose in the literature) due to the form of equation (19). Further details are supplied in the Appendix.
^{7}This would be possible for singlecrystal samples or polycrystalline samples with large coherently scattering domains (Saenger & Noyan, 2001).
^{8}The textbook by Smith is publicly available at https://www.analog.com/en/education/educationlibrary/scientist_engineers_guide.html. There are several typos in the text. Caveat emptor.
Acknowledgements
The authors thank Dr SeungYub Lee and Dr Shangmin Xiong for helpful discussions at the initial stages of this study.
Funding information
The following funding is acknowledged: Scientific and Technological Research Council of Turkey (TÜBİTAK) (award No. 118C268 to Hande Öztürk).
References
Avillez, R. R. de, Abrantes, F. G. & Letichevsky, S. (2018). Mater. Res. 21, e20170980. Web of Science CrossRef Google Scholar
Cline, J. P., Mendenhall, M. H., Ritter, J. J., Black, D., Henins, A., Bonevich, J. E., Whitfield, P. S. & Filliben, J. J. (2020). J. Res. Natl Inst. Stand. Technol. 125, 125020. Web of Science CrossRef Google Scholar
Dorofeev, G. A. A. N., Streletskii, I. V., Povstugar, A. V., Protasov, A. V. & Elsukov, E. P. (2012). Colloid J. 74, 675–685. Web of Science CrossRef CAS Google Scholar
Guinier, A. (1994). Xray Diffraction in Crystals, Imperfect Crystals, and Amorphous Bodies. Toronto: Dover Publications. Google Scholar
Ida, T., Shimazaki, S., Hibino, H. & Toraya, H. (2003). J. Appl. Cryst. 36, 1107–1115. Web of Science CrossRef CAS IUCr Journals Google Scholar
Kaszkur, Z. (2006). Z. Kristallogr. Suppl. 23, 147–154. CrossRef Google Scholar
Kaszkur, Z. (2019). J. Appl. Cryst. 52, 693–694. Web of Science CrossRef CAS IUCr Journals Google Scholar
Klug, H. P. & Alexander, L. E. (1974). Xray Diffraction Procedures for Polycrystalline and Amorphous Materials. New York: WileyInterscience. Google Scholar
Langford, J. I., Louër, D. & Scardi, P. (2000). J. Appl. Cryst. 33, 964–974. Web of Science CrossRef CAS IUCr Journals Google Scholar
Leoni, M. (2019). International Tables for Crystallography, Vol. H, edited by C. J. Gilmore, J. A. Kaduk & H. Schenk, ch. 5.1, pp. 524–537. Chichester: Wiley. Google Scholar
Lucks, I., Lamparter, P. & Mittemeijer, E. J. (2004). J. Appl. Cryst. 37, 300–311. Web of Science CrossRef CAS IUCr Journals Google Scholar
MathWorks (2019). MATLAB. Version 9.7R2019b. The MathWorks Inc., Natick, Masschusetts, USA. Google Scholar
Matyi, R. J. L. H., Schwartz, L. H. & Butt, J. B. (1987). Catal. Rev. 29, 41–99. CrossRef Web of Science Google Scholar
Mittemeijer, E. J. & Welzel, U. (2008). Z. Kristallogr. 223, 552–560. Web of Science CrossRef CAS Google Scholar
Noyan, I. C., Bunn, J. R., Tippett, M. K., Payzant, E. A., Clausen, B. & Brown, D. W. (2020). J. Appl. Cryst. 53, 494–511. Web of Science CrossRef CAS IUCr Journals Google Scholar
OriginLab (2021). OriginPro. Version 2021b. OriginLab, Northampton, Masschusetts, USA. Google Scholar
Palosz, B., Grzanka, E., Gierlotka, S. & Stelmakh, S. (2010). Z. Kristallogr. 225, 588–598. Web of Science CrossRef CAS Google Scholar
Patterson, A. L. (1939a). Phys. Rev. 56, 972–977. CrossRef CAS Google Scholar
Patterson, A. L. (1939b). Phys. Rev. 56, 978–982. CrossRef CAS Google Scholar
Saenger, K. L. & Noyan, I. C. (2001). J. Appl. Phys. 89, 3125–3131. Web of Science CrossRef CAS Google Scholar
Scardi, P. & Leoni, M. (2006). J. Appl. Cryst. 39, 24–31. Web of Science CrossRef CAS IUCr Journals Google Scholar
Schwartz, L. H. & Cohen, J. B. (1987). Diffraction from Materials. Heidelberg: Springer. Google Scholar
Smith, S. W. (1999). The Scientist and Engineer's Guide to Digital Signal Processing. San Diego: California Technical Publishing. Google Scholar
Taylor, A. (1961). Xray Metallography. New York: John Wiley and Sons. Google Scholar
Tomaszewski, P. E. (2013). Phase Transitions, 86, 260–266. Web of Science CrossRef CAS Google Scholar
Warren, B. E. (1969). Xray Diffraction. Reading: AddisonWesley. Google Scholar
Warren, B. E. & Averbach, B. L. (1950). J. Appl. Phys. 21, 595–599. CrossRef CAS Web of Science Google Scholar
Warren, B. E. & Averbach, B. L. (1952). J. Appl. Phys. 23, 497. CrossRef Web of Science Google Scholar
Williamson, G. K. & Hall, W. H. (1953). Acta Metall. 1, 22–31. CrossRef CAS Web of Science Google Scholar
Wolfram (2021). Mathematica. Version 12.3.1. Wolfram Research Inc., Champaign, Illinois, USA. Google Scholar
Xiong, S., Öztürk, H., Lee, S.Y., Mooney, P. M. & Noyan, I. C. (2018). J. Appl. Cryst. 51, 1102–1115. Web of Science CrossRef CAS IUCr Journals Google Scholar
Xiong, S., Öztürk, H., Lee, S.Y., Mooney, P. M. & Noyan, I. C. (2019). J. Appl. Cryst. 52, 695–696. Web of Science CrossRef CAS IUCr Journals Google Scholar
Ying, A. J., Murray, C. E. & Noyan, I. C. (2009). J. Appl. Cryst. 42, 401–410. Web of Science CrossRef CAS IUCr Journals Google Scholar
Young, R. A., Gerdes, R. J. & Wilson, A. J. C. (1967). Acta Cryst. 22, 155–162. CrossRef IUCr Journals Web of Science Google Scholar
This is an openaccess article distributed under the terms of the Creative Commons Attribution (CCBY) Licence, which permits unrestricted use, distribution, and reproduction in any medium, provided the original authors and source are cited.