X-ray powder diffraction in education. Part I. Bragg peak profiles

Peak profile functions in powder diffraction are presented with the support of Mathematica scripts, easily usable by interested readers, to explore the effect of different instrumental and microstructural parameters. Subsequent articles will illustrate the most common aberrations to the measured patterns and then total scattering methods.


Introduction
X-ray powder diffraction (XRPD), an established technique for the study of crystal structures, finds increasing application in qualitative and quantitative phase analysis and the study of microstructure for micro-and nanocrystalline materials.
For the study of crystalline materials, a powder pattern can be divided into Bragg reflections and background, where the latter among other contributions contains diffuse scattering from the sample. Bragg reflections are characterized by their position, intensity, breadth and shape, each containing a wealth of information (Dinnebier & Billinge, 2008;Dinnebier et al., 2018).
Computer programs for fitting peak profiles in powder diffraction perform either single peak fits or whole powder profile fitting with or without a structural model. Many programs have in common the ability to work like a black box for the user. It is frequently the case that some parameters are little affected by the nonlinear least-squares procedure, while others change wildly. In particular, among non-experts, a lot of confusion exists regarding the meaning and significance of fixed or refined parameters, their contribution to the Bragg peaks, the quality of the fit, the range of convergence, precision and accuracy, meaning and reliability of standard deviations etc. In addition, many functions describing a physical effect show discontinuities, and the refined parameters are typically restricted by physical boundaries. Understanding the algorithms is necessary for correct use of the software and assessment of the reliability of results.
In this series of papers dealing with the visualization of mathematical functions used to describe a powder pattern, we present a collection of user-friendly, interactive and freely distributable Mathematica (Wolfram Research, https:// mathworld.wolfram.com/) teaching scripts. All scripts have been written in Wolfram Mathematica version 12.1.1.0 and are constantly updated. They are freely available at the TOPAS Wiki web site (http://topas.dur.ac.uk/topaswiki). Non-subscribers of Mathematica can run the scripts using the freely available Wolfram Player at https://www.wolfram.com/player/. Bugs and problems should be reported to r.dinnebier@fkf.mpg.de. In particular, the so-called 'Manipulate' option of Mathematica is extensively used to visualize the impact of parameters in an interactive manner. When possible, parameters from real-life examples are used as the default. The idea is 'learn by doing', to gain intuition for what a mathematical model does to the diffraction peaks in a powder pattern and what the limitations of the said model are, since every model is an oversimplification of the underlying physics.
At the moment, the following three parts of this series are envisaged: part I, the peak profile of a powder pattern (this paper); part II, common correction functions in powder diffraction; part III, total scattering.
In this first paper, the various contributions to the peak profile in a powder diffraction pattern are described in terms of physical models and the effect they have on the peak shape. In general, the profile of a Bragg peak has contributions from the diffractometer and the microstructure of the sample. In the past few decades, a variety of sophisticated techniques and computer programs for analysing the peak shape of a powder pattern have been developed (Dinnebier & Billinge, 2008).
In general, the profile Èðx 0 Þ of a Bragg reflection centred at the peak position X 0 can be approximated by mathematical convolution (denoted by the symbol Þ of contributions from the instrument, called the instrumental resolution function (IRF), and from the microstructure (MS) of the sample (Klug & Alexander, 1974): where X is the observed peak position on the scale in which the data are recorded or analysed. The profile function is therefore described relative to the peak centre X 0 . Since the majority of powder patterns are directly measured as a function of the scattering angle 2 in degrees, we will use 2 as an independent variable for our profile functions throughout the text. Note that many of the aberrations to be discussed in part II are specifically 2 dependent. The conversion to reciprocal space in nm À1 is given by the Bragg equation: with the wavelength in nm and the distance between lattice planes d hkl with indices hkl in nm. Other scales besides 2, d and d Ã include Q (= 2d Ã ), time of flight (TOF) and energy (E). The conversion factor c ¼ 180 which follows from the differential of the Bragg equation, is used to change space from 2 in degrees to d Ã in nm À1 . [The factor /180 is used to convert an angular value given in degrees (e.g. the breadth of a reflection) to radians.] The Fourier transformation (FT) of a function in 2 space in degrees leads to a function with reciprocal angle and units of degree À1 that we will use for illustrative purposes only. The FT of a function in d Ã space and units of nm À1 leads to a function in nm on a length scale in real space. The IRF itself can be considered a convolution of contributions coming from the finite width of the X-ray source (X-ray tube or synchrotron), called the emission profile, and a series of horizontal and vertical instrumental aberrations due to the diffractometer. In the most popular configuration, these include the angular acceptance function of the Soller slit(s) controlling the axial beam divergence, the angular acceptance function of the plug-in slit controlling the equatorial beam divergence, the angular acceptance function of the receiving slit etc. For linear position-sensitive detector (PSD) systems, the receiving slit aberration is replaced by functions describing the defocusing due to asymmetric diffraction, the parallax error and the point spread function of the detector: The same principle holds for the MS contribution of the sample, which can be viewed as a convolution of contributions from effects like the size of coherently scattering domains, isotropic and/or anisotropic microstrain, faulting etc.: Equations (5) and (6) implicitly assume that all components of the convolution are independent, which is a common approximation but not exact for MS parameters.
To describe the peak shape in a powder diffraction pattern, two approaches are common nowadays: (i) Ideally, all contributions and their mathematical description are known. The resulting peak shape can thus be calculated from first principles. For building the IRF, the fundamental parameters (FP) approach (Cheary & Coelho, 1992, 1998 and for the MS the whole powder pattern modelling (WPPM) concept (Scardi & Leoni, 2001 are quite common. (ii) Alternatively, because many contributions to powder diffraction peaks have a nearly Gaussian or Lorentzian shape, the Voigt function, which is a convolution of Gaussian and Lorentzian components, or more commonly the pseudo-Voigt function for faster computation, is widely used to describe peak profiles in powder diffraction. Importantly, the latter approach is empirical, and all interpretations of MS parameters based on the Voigt profile must be carefully evaluated.
In fact, real instruments are always more complex than any possible FP model. The differences between model and reality can be a source of systematic errors with the risk of ignoring some features of the IRF. As a commonly accepted compromise, the IRF is obtained by fitting a line profile standard with virtually no MS contributions using an empirical Voigt profile, eventually convoluted with some instrumental aberration functions. In contrast, the MS is derived from more physical or teaching and education at least phenomenological functions with few but meaningful parameters.

Fourier transformation and convolution
The FT greatly simplifies the convolution process and reveals many important properties of the peak profile, making it a good mathematical procedure to begin with.
The FT of a function f ðxÞ is a complex function defined as The Fourier integral transforms a function f ðxÞ (x in real space) by an integral over cosine and sine functions to FðsÞ (s in reciprocal space = Fourier space). Some common choices for the Fourier parameters {a, b} which have a big impact on scaling are {0, 1} (modern physics), {1, À1} (pure mathematics, systems engineering), {À1, 1} (classical physics), {0, À2} (signal processing) and {0, 2} (crystallography).
The back-transformation is given by The independent variables, s and x, are reciprocal to one another, so the exponent of the exponential function is dimensionless. If x is in position space, s is in 'Fourier space' or 'reciprocal space'. In the following, we will use the crystallographic Fourier parameters {0, 2} for simplicity: Convolution, or folding, is a basic concept in crystallography that is particularly important in powder diffraction analysis. The process of convolution is one in which the product of two functions f ðxÞ and gðxÞ is integrated over all space: where hðx 0 Þ is the convolution product, y is the variable of integration in the same domain as x 0 and denotes the convolution process. Convolution can be understood as 'blending' one function with another, producing a kind of very general 'moving average' [see Weisstein (2021) for a definition and animated examples]. Most functions cannot be convoluted analytically and the convolution integral needs to be calculated numerically.
An alternative method of calculation follows from the convolution theorem of FT. The FT of the convolution function can be calculated by which can be rewritten using the substitution u ¼ x 0 À y and therefore du ¼ dx 0 as f ðyÞgðuÞ expð2isyÞ dy expð2isuÞ du f ðyÞ expð2isyÞ dy gðuÞ expð2isuÞ du From the last, it follows directly that the FT of the convolution integral is the product of the FTs of all functions participating in the convolution, while the back-transformation of a convolution is the product of the back-transformed functions which participate in the convolution: In many cases, an analytic FT of a mathematical function is not possible, requiring the numerical FT of a finite list of N function values uðrÞ at index r (running from 1 to N), which are separated by a constant step width Áx. The corresponding running index in Fourier space is s (running from 1 to N) and the FT values are separated by a constant step width. This means that some piece of the integral can be discretized by x ! x 1 þ ðr À 1ÞÁx and s ! ðs À 1ÞÁs In Mathematica, a numerical discrete FT is defined as F s ¼ ð1=N 1=2 Þ P N r¼1 u r exp½2iðr À 1Þðs À 1Þ=nÁx, from which it follows that the f r values must be multiplied by a complex phase factor u r = f r N 1=2 exp½2iðs À 1Þx 1 =ðNÁxÞ (Mathematica stackexchange questions 1714, https://mathematica. stackexchange.com/questions/1714/numerical-fourier-transformof-a-complicated-function/151179#151179).

Probability density functions
The shape of powder diffraction peaks is usually described by individual or convoluted probability density functions (PDFs) used to describe statistical processes. In statistics, PDFs are characterized by their moments, which are quantitative measures related to their shape. The zeroth moment is the total probability (which equals one), the first moment is the expected value (mean), the second central moment is the variance, the third standardized moment is the skewness and the fourth standardized moment is the kurtosis. The variance is the expectation of the squared deviation of a random variable from its mean. Skewness is a measure of the asymmetry of the probability distribution of a real-valued random variable about its mean, and kurtosis is a measure of the 'tailedness' of the probability distribution of a real-valued random variable (from Wikipedia).
In powder diffraction, the widths of distributions used to describe peak profiles are either the standard deviation (), the full width at half-maximum (FWHM) which is the difference between the two values of the independent variable at which the dependent variable is equal to half of its maximum value (height of the peak), the half-width at half-maximum (HWHM) which is FWHM/2, or the integral breadth (), the last being the width of a rectangle of the same height and area as the peak.
The standard deviation or square root of the variance of a distribution described by f(x) (which is normalized to an area of unity) is given as and the average is given as which agrees (in the case of a symmetric unimodal distribution) with the peak maximum x 0 . For many common distributions, the variance and therefore the standard deviation are not defined. In these cases, the HWHM, FWHM or is always used. The last is defined as for a function f(x) that is normalized to an area of unity.

Gaussian distribution
The Gaussian (or normal) distribution is a very common distribution. Physical quantities expected to be the sum of many independent processes (such as measurement errors) often have distributions that are nearly normal. A Gaussian function normalized to an area of unity centred at a mean or expectation value X 0 is ( Fig. 1, left) where is the standard deviation and 2 the variance. The pre-factor follows from normalization: The HWHM of a Gaussian is given by With HWHM G as the width parameter, the definition of the Gaussian changes accordingly to The integral breadth of a Gaussian is given by The normalized FT (maximum value is 1) of a Gaussian function is itself a Gaussian function (Fig. 1, right): with an FWHM of The convolution of a Gaussian (1) with another Gaussian (2) is a Gaussian with the following property: which follows directly from the additivity of variances and also holds for integral breadths.

Lorentz distribution
The Lorentz (or Cauchy) function is another important continuous probability density distribution, which might be attributed, for instance, to the lifetime broadening of the characteristic X-ray emission line. The Lorentz function is used, for example, to describe the emission profile from an X-ray tube and the peak profile in some faulting problems when a random probability exists in the stacking sequence, as well as to approximate crystallite size and strain effects from the sample. Also, in a perfect infinite crystal, Bragg peaks are not functions but finite Lorentzians with the FWHM being the Darwin width (Warren, 1990). The Lorentz distribution normalized to unity is defined as (Fig. 2, left) with being the Lorentzian HWHM HWHM L . The Lorentz distribution is an example of a distribution with no mean, no standard deviation, and no variance or higher moments defined. Its mode and median (the value separating the higher half from the lower half of a data sample) are well defined and are both equal to X 0 . The integral breadth of a Lorentzian is given as The real part of the FT (normalized to unity) of a Lorentzian is ( Fig. 2, right) with an FWHM of The convolution of a Lorentzian (1) with another Lorentzian (2) is a Lorentzian with the following property: which also holds true for integral breadths.

The Voigt distribution
The Voigt distribution, named after the German physicist Woldemar Voigt (Voigt, 1912), can be regarded as the convolution of a normalized Gaussian and a normalized Lorentzian: No analytical solution exists for the convolution integral, but it can be expressed by the real part of the complex error function for which good approximations exist (Fig. 3, left): where wðzÞ is called the Faddeeva function (also called the Kramp function or relativistic plasma dispersion function) and is a scalable complex conjugated error function. This is given by teaching and education  wðzÞ ¼ expðÀz 2 Þ erfcðÀizÞ and the argument z is When the results of a series of measurements are described by a normal distribution with standard deviation and expected value 0, then is the probability that the error of a single measurement lies between Àa and +a, for positive a. The error function is defined as (Fig. 4) The complementary error function is defined as The FWHM V of a Voigt function can be reasonably approximated either by (Olivero & Longbothum, 1977) or by a polynomial of fifth order (Thompson et al., 1987): The integral breadth of a Voigt function is given as (Schoening, 1965) with the characteristic breadth ratio of a Voigtian k given by The FT of the Voigt function can be written as (Fig. 3,   Series of normalized Voigt profiles (left) and their corresponding real Fourier transforms (right) with different Gaussian standard deviation and Lorentzian HWHM on a 2 and reciprocal 2 scale.

Figure 4
Plot of the error function and the complementary error function given in equations (38) and (39) [from Dinnebier et al. (2018)].

Figure 5
Plot of the FWHM of the TCH pseudo-Voigt function as a function of the FWHMs of the Gaussian and the Lorentzian part (solid lines) and the corresponding integral breadths (dashed lines).
with an FWHM of The convolution of two Voigt functions again is a Voigt function, a property that is used, for example, in the double-Voigt approach where two Voigtians are convoluted, each having its own Gaussian and Lorentzian fractions (Balzar, 1999). For a long time, the exact computation of a Voigt profile was computationally expensive. One way to approximate a given Voigt function is to use a linear combination of a Gaussian and a Lorentzian, called the pseudo-Voigt approximation: The mixing parameter can be calculated using equation (43) in the range 0 1 by a cubic polynomial leading to the Thompson-Cox-Hastings (TCH) pseudo-Voigt function (Thompson et al., 1987): If a purely empirical approximation of a Voigt function is needed (as is generally sufficient for the determination of the IRF), a simplified version of the pseudo-Voigt function can be used: with its FT Among the many approximations of the Voigt function found in the literature, the one by Chiarella & Reichel (1968 Screenshot of a Mathematica script for comparing the pure Voigt profile, the Reichel approximation and the TCH pseudo-Voigt profile for a given Gaussian standard deviation and Lorentzian HWHM in 2 or d Ã space. Values for crystallite size and microstrain were calculated using the formulas by Balzar (1999) under the assumption that the entire width of the peak is caused by either one. A good explanation of how size and strain values are calculated from the FWHM or the integral breadth of a Voigt function in the TOPAS (Coelho, 2018) software is given by Evans (2021).
which was later modified by Abrarov et al. (2012), is a good compromise regarding speed and accuracy (Figs. 6 and 7): with h ¼ = m , and the chosen parameters m and N, which determine the degree of approximation, typically m = 12 and N = 23. This function can be efficiently compiled in Mathematica (Mathematica stackexchange question 19884, https:// mathematica.stackexchange.com/questions/19884/compilingthe-voigtdistribution-pdf), increasing the speed of computation considerably. Since the difference between a Voigt profile and the Chiarella and Reichel approximation is negligible (Fig. 7), the latter is used throughout the Mathematica scripts if a Voigt profile is required. The dependence of the Gaussian FWHM G and Lorentzian FWHM L on the 2 diffraction angle is usually described by low-order polynomials as a function of cos À1 and tan . While the cos À1 dependence follows directly from the Scherrer equation (Scherrer, 1918) and is attributed to the size of coherently scattering domains (often denoted as crystallite size), the tan dependence is a measure of microstrain and can be derived from the full derivative of the Bragg equation (Bragg & Bragg, 1913). In an empirical but flexible way, the Gaussian FWHM G and Lorentzian FWHM L are often defined as (Thompson et al., 1987;Young, 1993) where U, V, W, X, Y, Z are refinable parameters. Apparently, the parameters have been chosen in such a way that U and X are somehow related to microstrain while Z and Y are kind of related to domain size (Von Dreele, 2008) (Fig. 8). Refining all parameters simultaneously generally leads to overparameterization and high correlation, destabilizing the optimization process.
In order to retrieve the width parameters of a Voigt function which is blurred by Poisson noise, a variety of leastsquares algorithms (e.g. gradient, Levenberg-Marquardt, Newton) are available (Press et al., 2007). A Mathematica script has been written to evaluate the reliability of parameter retrieval depending on the amount of noise and to determine the degree of correlation between the fitting parameters ( Fig. 9). The script is a useful tool to test and visualize data quality and the effect of noise on the refined parameters, applying different algorithms of minimization. It also introduces the concept of estimated standard deviation.

Instrumental aberrations to the peak profile
Only in a few cases can the IRF be fully described using a symmetrical Voigt-like profile function. In particular, the Screenshot of a Mathematica script for comparing the Reichel approximation of the Voigt function with the TCH pseudo-Voigt profile as a function of Gaussian standard deviation and Lorentzian HWHM at a given position relative to the peak centre in three dimensions.

Figure 8
Screenshot of a Mathematica script for visualization of the dependence of the TCH pseudo-Voigt function and its Fourier transform on the diffraction angle. Values for crystallite size and microstrain were calculated using the formulas by Balzar (1999) under the assumption that the entire width of the peak is caused by either one.

Figure 9
Screenshot of a Mathematica script for visualization of least-squares fitting using different algorithms of the TCH pseudo-Voigt perturbed by Poisson noise.
curvature of the Debye-Scherrer rings introduces a certain amount of asymmetry due to axial divergence if a rectangular receiving slit or silicon strip is cutting out only a small part of the ring. Each element in the optical pathway of the incident and diffracted beam convolutes its characteristic shape into the peak profile. The functions to describe the effects are either highly sophisticated based on physics and geometry, purely empirical, or a mixture of both. They are convoluted in real space, or multiplied in Fourier space, into the existing peak profile. In the following, a selection of four simple aberration functions in powder diffraction, which have proven to be useful, are discussed. Although this section is restricted to instrumental aberrations, the asymmetry introduced by the transparency effect, which depends on the preparation technique, is included since it is not related to the MS of the sample. Note that all functions presented here can be found in the TOPAS software (Coelho, 2018).

The box (hat) function
Several instrumental aberrations in the equatorial plane of a diffractometer are commonly described by a box (sometimes called top-hat) function. These include the width of the source, the thickness of the sample surface as projected onto the equatorial plane, the width of the receiving slit, the width of strips in position-sensitive strip detectors etc. The normalized (area under its graph is 1) box function with width a is defined as (Fig. 10, left In the extreme case of a ! 0, the box function turns into a function. For practical reasons this is achieved by setting a to a value <10 À5 . The normalized (maximum value is 1) FT of a box function with the reciprocal variable s is a real function calculated as (Fig. 10, right) The convolution of a box function with itself leads to the triangular function. If the convolution process is repeated, the convoluted function approaches a Gaussian. In practice, this is realized after five or more convolutions.
To mimic the transmittance of a rectangular slit or a source with a width of c in mm (typically on the order of 0.1 mm) in the equatorial plane of a diffractometer with a secondary radius R s in mm, a constant function for calculating the FWHM a of the box function is used (Fig. 11): In order to calculate the width (FWHM) of the box function in degrees for the specimen tilt, the following function is used (Fig. 11): where c (in mm) represents the 'thickness' of the sample surface as projected onto the equatorial plane (Cheary & Coelho, 1992). In the case that an equatorial aberration is non-symmetrical (as for tube tails which might have different lengths on each side), a normalized halfbox function acting either on the left or the right side of a diffraction peak can be defined as (Fig. 12  The box function on a 2 scale with a width of 0.5 2 (left) and its real Fourier transform (right).

Figure 11
Aberrations in the equatorial plane of a diffractometer, assuming a secondary radius of 217.5 mm and a slit width (either receiving slit or size of the source) and a specimen tilt of 0.1 mm, as a function of diffraction angle 2.
with its normalized complex FT (Fig. 12, right) To mimic the aberration caused by tube tails (Bergmann et al., 2000), a combination of a box function describing the tube filament width in mm and two halfbox functions with an effective width of tube tails in the equatorial plane perpendicular to the X-ray beam in the negative and positive z direction in mm can be used. An additional parameter defining the fractional height of the tube tails relative to the main beam is then needed (Fig. 13). A screenshot of a Mathematica script where the box/halfbox function from different kinds of aberrations in the equatorial plane of the diffractometer is convoluted into a Voigt profile is shown in Fig. 14.

The circles function
A simple approximate function for modelling the asymmetry of a Bragg reflection is the circles function, with the cut-off value " m determining the curvature as an adjustable parameter (Fig. 15, left): " m for 0 < x 0 " m and j" m j > 0: One of the main applications of this function is the phenomenological modelling of the peak asymmetry caused by axial divergence. This is predominantly due to the increasing curvature of the Debye-Scherrer rings at very low and extremely high angles which are cut by (typically) rectangular receiving slits of finite width (Cheary & Coelho, 1998). The complex FT of the circles function (Fig. 15, right) is defined as CIRCLESðsÞ ¼ 2½ 1 F 2 ð0:25; 0:5; 1:25; À0:25s 2 " 2 m Þ þ 0:3325is" m 1 F 2 ð0:75; 1:5; 1:75; À0: with the hypergeometricPFQ function 1 F 2 . In order to model axial divergence with the circles function for angular dispersive data, the angular dependence of the cutoff value " m is well described by a tan (2) À1 type of function of the receiving slit length c in mm. Note that the asymmetry is reversed above 90 2 [actually closer to 120 2 as described by Cheary & Coelho (1998)] (Fig. 16): More sophisticated functions to describe the asymmetry due to axial divergence of divergent beam diffractometers can be found in the literature. The mathematical formalism to describe peak asymmetry due to the finite size of the detector receiving slit and the curvature of the Debye-Scherrer rings for parallel beam geometry has been developed by van Laar & Yelon (1984) and was implemented by Finger et al. (1994). Extensions to cover divergent beam geometry, common for laboratory diffractometers, have also been developed (e.g. Cheary & Coelho, 1998;Mendenhall et al., 2017). The halfbox function on a 2 scale with a width of 0.25 2 (left) and its complex Fourier transform (right).

Figure 13
Tube tail function as a combination of a box function and two (left, right) halfbox functions for the tails.
A screenshot of a Mathematica script where the circles function for describing asymmetry due to axial divergence of the diffractometer is convoluted into a Voigt function is shown in Fig. 17.

1/x function
Another aberration function to describe certain types of asymmetry of a Bragg reflection is the 1/x decay function, defined as (Fig. 18, left) oneonx ðx 0 Þ ¼ 1 where the parameter " m , which can be either positive or negative, is the cut-off value defining the extension of the asymmetric tail on the given length scale. The normalized complex FT of a 1/x function with the reciprocal variable s is calculated as (Fig. 18, right) ONEONX ðsÞ ¼ 1:25 The circles function on a 2 scale with three cut-off values of 2, 4 and 8 2 (left) and the Fourier transform of the circles function with a cut-off value of 4 2 showing the real and the imaginary part (right).

Figure 14
Screenshot of a Mathematica script dealing with the convolution of the box function into a Voigt profile. Different kinds of aberrations in the equatorial plane of a diffractometer assuming a secondary radius of 217.5 mm (originating from a receiving slit, the size of the source or the specimen tilt) as a function of diffraction angle 2 can be selected.
with the two Fresnel integrals The 1/x function is mainly used for the phenomenological description of the equatorial divergence caused by a divergence slit (either fixed or variable) on the peak profile for angular dispersive data. The cut-off value " m as a measure for the peak asymmetry is calculated for a fixed or variable divergence slit as (Cheary et al., 2004)  Cut-off values (= curvature) of the circles function to describe asymmetry due to axial divergence as a function of scattering angle for a diffractometer with a secondary radius of 217.5 mm and different receiving slit lengths in mm.

Figure 17
Screenshot of a Mathematica script dealing with the convolution of the circles function for modelling asymmetry due to axial divergence caused by the length of a rectangular receiving slit.

Figure 18
The 1/x function on a 2 scale with three cut-off values of 2, 4 and 8 2 (left) and its Fourier transform with a cut-off value of 4 2 showing the real and the imaginary part (right).
opening angle (in degrees) if the divergence slit is fixed (Fig. 19): A screenshot of a Mathematica script where the 1/x function for describing asymmetry due to fixed or variable divergence slits of the diffractometer is convoluted into a Voigt function is shown in Fig. 20.

The exponential function
Asymmetry of Bragg reflections can also be described by a normalized exponential asymmetry decay function of the type (Fig. 21, left) [ln(0.001) is a scaling factor which ensures that the value at x 0 ¼ " m is 0:001=0:14462j" m j] that is convoluted into the peak profile. The parameter " m , which can be either positive or negative, is the cut-off value which defines the extension of the asymmetric tail on the given length scale and is a measure for the degree of asymmetry. The normalized complex FT of the expo function with the reciprocal variable s is calculated as (Fig. 21, right) EXPO ðsÞ ¼ ð2Þ 1=2 ð2:7586iÞ À ð0:0027586iÞ expði" m sÞ À ln 0:001i þ " m s : The exponential function can be useful for describing, for example, the highly asymmetric instrumental peak shape of TOF data or the effects of transparency on the peak shape in Bragg-Brentano geometry.
The main application of the exponential aberration function is to describe the peak asymmetry caused by the transparency effect, where a low-absorbing sample is filled in a deep cavity Screenshot of a Mathematica script dealing with the convolution of the 1/x correction function for modelling the peak asymmetry due to equatorial divergence caused by fixed or variable divergence slits.

Figure 19
Angular dependence of the 1/x correction function for modelling the peak asymmetry due to equatorial divergence caused by a fixed divergence slit of 0.5, 1.0 or 2.0 opening (left) and a variable divergence slit with a constant illuminated sample length of 10, 20 or 30 mm (right). The secondary diffractometer radius is set to 217.5 mm.
of a flat-plate sample holder in Bragg-Brentano geometry. The diffracted peak thus is a convolution from multi-diffraction at different depths with increasing absorption.
The cut-off value " m as a measure for the peak asymmetry is calculated for a linear absorption coefficient a in cm À1 (Fig. 22) as ð2Þ ¼ 900 ln 0:001 sin 2 aR s : A screenshot of a Mathematica script where the exponential function for describing asymmetry due to transparency is convoluted into a Voigt function is shown in Fig. 23.

Microstructural contribution to line profiles
Once determined, the instrumental contribution to the diffraction profile can be used constantly for a given instrument and setup. Whether it is modelled on the basis of The exponential function on a 2 scale with three cut-off values of 2, 4 and 8 2 (left) and its Fourier transform with a cut-off value of 4 2 showing the real and the imaginary part (right).

Figure 22
Angular dependence of the exponential function for modelling the peak asymmetry due to absorption for different absorption coefficients of 10, 50 and 100 cm À1 . The secondary diffractometer radius is set to 217.5 mm.

Figure 23
Screenshot of a Mathematica script dealing with the convolution of the exponential correction function for modelling the peak asymmetry due to the transparency effect caused by low linear absorption in flat-plate sample holders with sample-filled deep cavities.
fundamental parameters or with empirical approaches, the instrumental profile can be acquired experimentally, using appropriate powder standards (Cline et al., 2013). Separately, the MS contribution varies for each case study and, therefore, there cannot be a universally valid model. This explains the practical use, described before, of empirical profile functions with adaptive parameters, justified by the fact that they 'work well', i.e. they fit the patterns in most cases. However, the interpretation of these parameters can be misleading due to the substantial arbitrariness in the choice of profile functions and their combination (Scardi et al., 2004;Scardi, 2020).
The WPPM provides a description of the MS component of the profile based on a convolution of terms referring to physical parameters. The two normally prevailing effects -size and strain broadening -are introduced below, whereas other effects and more extensive discussions are reported in the literature (Scardi et al., 2018;Scardi, 2020, and references therein). This choice corresponds to writing equation (6) as where the size is expressed as a function of lattice direction and strain as a function of both lattice direction and d Ã . This highlights an important property: when the diffraction profile is represented in reciprocal space, the size component, unlike the strain component, is independent of d Ã , which allows for separation of the two effects, e.g. through the classic Warren & Averbach method (Warren & Averbach, 1950;Warren, 1990). Both components, however, can depend on Miller indices (hkl), respectively, for the shape of the crystalline domains, and for the anisotropy of defect strain field and of the elastic medium. Note that, for small nanoparticles/domains, the separation of domain size and strain parameters becomes increasingly difficult, and it becomes ill-defined for very small or highly disordered domains in a crystallographic sense. Within the limits of the tangent plane approximation (TPA) (von Laue, 1926;Beyerlein et al., 2011), it can be shown that the FT of the profile component of the size effect, A S (s), has a simple geometric interpretation (Wilson, 1962) (Fig. 24). In fact, it corresponds to the volume of intersection V(s, hkl) between the crystalline domain and the same domain translated by a distance s along [hkl], normalized to the domain volume, V 0 , Simple geometrical shapes allow for an easy evaluation of V(s, hkl) in closed analytical form (Scardi & Leoni, 2001). Such is the case for frequently observed nanocrystal shapes, like the cube, tetrahedron, octahedron and sphere. The last is frequently assumed, even if only approximately valid, when nanocrystals do not have a specific shape but are reasonably equiaxed. Then the FT of the line profile given by spherical domains of diameter D is (Fig. 25,    Geometrical interpretation of the Fourier transform of the size-effect profile component. The common volume function, V(s, hkl), along two different directions, [h 1 k 1 l 1 ] and [h 2 k 2 l 2 ] (left), with corresponding Fourier transforms (right); s max stands for the maximum dimension in the crystalline domain along the given direction.
Apart from the spherical shape, A S is generally a function of hkl, a dependence requiring the choice of lattice orientation within the crystalline domain. Even if complex shapes and low-symmetry structures require numerical calculation methods, the calculation of A S can lead back to a tractable geometry problem (Leonardi et al., 2012) (Fig. 25, right). It is unlikely that the domains, despite having the same shape, are of the same size. The problem can be addressed by considering a size distribution, as discussed by Scardi & Leoni (2001). A variety of distribution functions can be employed, but a lognormal distribution is often an appropriate choice (Fig. 26): The lognormal mean () and lognormal variance ( 2 ) give the central moments, M l;n , where from which the mean and standard deviation, respectively, are This distribution refers to a single dimensional parameter, D, the diameter of the sphere or edge of a cube, tetrahedron, octahedron etc. It is theoretically possible to use more distributions, for example, base diameter and height for a cylindrical shape, but the strong correlation between the different parameters (two for each distribution) makes the use impractical and refinement by least-squares usually unstable. For a lognormal distribution of spheres, the FT A S l ðsÞ is (Scardi et al., 2018) A S l ðsÞ ¼ q 0 ðsÞ À 3 2 q 1 ðsÞ þ 1 2 q 3 ðsÞ ð 78Þ with q n ðsÞ ¼ s n erfc ln s À À ð3 À nÞ 2 where the H c n and K c coefficients depend on the domain shape (e.g. H c 0 = 1, H c 1 = À3/2, H c 2 = 0, H c 3 = 1/2, K c = 1, for spheres), and in particular K c ðhklÞ sets the maximum dimension in the crystalline domain along the given [hkl] (the formulas for the cubic crystal system require the condition h ! k ! l) [see Leonardi et al. (2012) for other shapes] (Fig. 27).
The effect of a lognormal size distribution of particles with a given shape and size on the pseudo-Voigt profile and its FT is shown in Fig. 28.
The WPPM approach for the size effect is flexible enough and appropriate to deal with many real cases. It generally Lognormal distribution of sizes with lognormal mean = 1.95 and lognormal standard deviation = 0.35, leading to a mean size of 7.47 nm (red arrow) and a standard deviation of s.d. = 2.68 nm.

Figure 27
Fourier transforms of the line profile of a lognormal system of spherical (1), cubic (2), tetrahedral (3) and octahedral (3)  works if domains are not too small, because the TPA does not hold for very small sizes ($5 nm is a realistic limit), and do not involve 'non-crystallographic' shapes like decahedra or icosahedra, which break translational symmetry. The strain effect is more complex and varied than the size effect, and can be treated in a simple way only as a perturbation of otherwise perfect crystalline domains. In these terms, it can be shown that the FT of the strain component, A D , of the line profile can be approximated by (Warren & Averbach, 1952) A D ðs; hklÞ ffi exp½À2 2 d Ã2 hkl hðÁsÞ 2 hkl i; where hðÁsÞ 2 hkl i is the variance of the displacement distribution for any couple of scattering centres in the domain at a distance s along the direction [hkl]. This can also be written in terms of 'microstrain', defined as h 2 hkl i 1=2 = hðÁsÞ 2 hkl i 1=2 =s. Expressions for the microstrain can be derived for specific case studies, e.g. the presence of dislocations (Wilkens, 1970a,b), or by adopting phenomenological models (Adler & Houska, 1979). In both cases, basic features are (i) a marked dependence on the reciprocal-space vector [d Ã2 hkl in equation (68)] and (ii) an anisotropy specific to the elastic properties of the material and of the type of lattice defects.
If one follows the Krivoglaz-Wilkens (KW) (Wilkens, 1970a,b) approach for dislocation strain broadening, the following expression results: with as the average dislocation density, b the Burgers vector modulus for the given slip system and R e the effective outer cut-off radius of the dislocation strain field; f* is a known (Wilkens) function of s/R e (Wilkens, 1970a), and À hkl depends on the direction, in terms of a fourth-order invariant expres-sion of the Miller indices. In the most general case (triclinic), this reads (Scardi et al., 2018) where a is the first unit-cell parameter according to the crystallographic conventions. Coefficients E 1 , E 2 , . . . , E 15 can be calculated for specific strain fields, like that of dislocations, given the slip system and elastic constants (Martinez-Garcia et al., 2009). The anisotropy factor, À hkl , for dislocation strain is also referred to as the average dislocation (or orientation) factor, " C C hkl . Symmetry reduces the number of terms in equation (83) (Scardi et al., 2018) down to two for the most symmetric, cubic case, for which equation (83) simplifies to Values of A (= E 1 ) and B [= 2(E 2 À E 1 )] can be calculated (Martinez-Garcia et al., 2009) but are also known in parametric form for a wide range of anisotropy (Zener) ratios and dislocation types in cubic systems (Ungá r et al., 1999;Dragomir & Ungá r, 2002).
As already pointed out, it is also possible to use empirical expressions, such as that named after Popa-Adler-Houska (PAH), which extends the Adler & Houska (1979)  Screenshot of a Mathematica script dealing with the convolution of the size effect into a Voigt profile and its Fourier transform for a given lognormal distribution of a chosen shape.
Formally, the anisotropy factor is the same as in equations (82) and (83), but in this case coefficients E 1 , E 2 , . . . , E 15 become just fitting parameters, to be optimized together with and . Fig. 29 shows examples of applications of both models for two typical case studies, nanocrystalline body-centred-cubic (b.c.c.) iron and face-centred-cubic fluorite, respectively. It is interesting to note in these two examples, beyond the use of the two different models of equations (82) and (85), the opposite trend of the anisotropy factor. B in equation (84) is negative for b.c.c. iron (like most metals) and positive for fluorite (like many binary salts), which means that the stiff/soft directions are opposite: [hhh]/[h00] for iron and [h00]/[hhh] for fluorite. As a consequence, (hhh) line profiles in iron tend to be less broadened than the (h00) profiles, apart from the d Ã hkl dependence, and vice versa for fluorite. Strain broadening anisotropy definitely helps the separation of size and strain effects.
MS effects can be visualized by two common representations, the Williamson-Hall plot (Williamson & Hall, 1953) and the Warren plot (Warren & Averbach, 1950). The former, consisting of a plot of the integral breadth (IB = peak area/ peak maximum intensity) versus d Ã hkl for as many peak profiles as are available in the experimental pattern, is also used for a preliminary assessment of the size and strain parameters. Note, however, that this analysis (also known as the Williamson-Hall method) is affected by the arbitrary choice of additivity for the size and strain IBs (Scardi et al., 2004). On the basis of the WPPM approach, instead, we can calculate the IBs from the FTs, for individual profile components and their convolution: where the ellipsis indicates that other terms contributing to the profile can be added (e.g. the instrumental profile). The dependence on hkl is invariably present in the strain term, whereas it is absent in the size effect for the frequent case where a spherical domain shape (average) is adopted. IBs of individual profile components are also easily obtained (in analytical form in some cases) as Analogous expressions hold for other possible profile components. All integrations virtually extend from À1 to +1, but more realistically within a finite range, determined by the maximum length for a given domain shape and the dispersion of the size distribution. Fig. 30 shows examples of the Williamson-Hall plot where the different dependence on d Ã hkl and on hkl can be seen. The different trends, with the constant size term and the characteristic anisotropy of the strain term, further demonstrate the possibility of separating the different effects.  The Warren plot concerns only the strain component of the profile and provides a convenient visualization of the dependence on s, the Fourier variable (a physical length in the crystalline domain), and anisotropy. It is a plot of ½hðÁsÞ 2 i hkl 1=2 , the standard deviation of the displacement distribution, as a function of s, for different directions [hkl]. Examples are shown in Fig. 30.