Using the BWA (Bertaut-Warren-Averbach) Method to Optimize the Technology of Crystalline Powders such as LiFePO4

1. Introduction

P. Scherrer began studying the shape of X-ray lines (LPA - line-profile analysis) to determine the sizes of nanometer objects in 1918. The history is described in [1, 2], and the latest results in [3]. In 1950, Bertaut [4] and Warren, Averbach [5] showed that the total diffraction intensity is the sum of individual diffractions from columns perpendicular to the planes that make up the crystallite. The magnitude of the size-induced broadening of diffraction reflections depends not only on the average sizes of crystallites, but also on the statistics of their distributions. The situation is described in detail in [4-9] and the following equation is assumed to hold

$$D_{i XRD}\equiv {\overline{D}}_{i }\left(M_i\right) = \frac{1}{V}\int M_i dV= K_n{\overline{L}}_{i } ,$$

(1)

where $M_i$ is the length of the column along the i-th direction, V is the crystallite volume, $K_n$ is the normalization coefficient of value $D_{i XRD}$ per crystallite linear size along the i-th axis. The result of averaging for a single ellipsoid, i.e. for a powder with a small dispersion, is obtained by integration [6-9].

Thus, the task of using the BWA method is reduced to determining the normalized coefficient $K_n$ by integrating (1) and summing over all crystallites, taking into account their shape, as well as the statistics of their distributions.

Based on TEM and SEM studies [2, 7, 10, 11], various size distribution functions: Lognormal, Normal, Gamma, and Poisson are normally used, more precisely, their one-dimensional projections onto some linear dimension. The first two of them can be used in their 3-dimensional versions with pairwise correlations between sizes along the coordinate axes [12-14], while for the last two such a possibility is not “transparent and easily accessible” [15-17]. In some exceptional cases, analytical approximations of distributions with a large dispersion of sizes may be useful [11]. This is probably also due to the fact that in some ranges of dispersion and average sizes, the Lognormal distribution is visually close, for example, to the Gamma distribution (see SI).

Currently, a large number of programs for processing diffractograms and determining $D_{i XRD}$ have been developed. First of all, it is necessary to take into account the instrumental contribution to diffraction line width [18, 19], and then the general microstructure refinement of all peak profiles (WPPM – Whole Powder Pattern Modeling) procedure can be used [3, 10, 11]. The size–strain line broadening analysis can be combined with evaluating the crystallographic texture of various nanosized powders [20,21,22] in the Rietveld refinement software known as Material Analysis Using Diffraction (MAUD) [23]. In TOPAS (Bruker) software anisotropic crystallite size analysis is also implemented. The accuracy of the double-Voigt approach and the dimensional parameters obtained for triaxial ellipsoid, elliptic cylinder and cuboid is discussed in [24, 25].

Thus, there is a fairly reliable procedure for determining the X-ray diffraction sizes of anisotropic powders, in which their size distribution is assumed to be known a priori, or by inspecting histograms obtained by counting a small number of particles. In the absence of more complete information about the 3-dimensionality of these distributions, the use of the BWA method in the form of equation (1) loses its meaning due to a large error. The known value of the coefficient $K_n$ can reduce the error that is very useful for more precise control of technology and predicting the target parameters of powders.

In this work, this problem is solved by the following numerical calculations: 1) the histograms of particle size distributions are simulated using 3D Lognormal or Normal functions in the form of a 3-dimensional and N-bit matrix, 2) the matrix of averaged columns of crystallites is calculated with the shape of g-type superellipsoids and the sizes corresponding to the histogram bins, 3) the value of ${\overline{D}}_{i }\left(M_i\right)$, defined as the average column length of the whole powder sample along the i-th axis, is calculated by element-by-element multiplication of these two matrices and then summing the matrix elements along the i-th axis.

2. Calculation model

2.1. Shape of crystallites

Lamé's surface of g-type superellipsoids is described by the implicit equation

$${\left(\frac{2 x_3}{L_3}\right)}^g+{\left(\frac{2 x_2}{L_2}\right)}^g+{\left(\frac{2 x_1}{L_1}\right)}^g=1,$$

(2)

superelipsoids belong to the superquadric family [26]

$${\left({\left|\frac{2 x}{L_x}\right|}^{\frac{2}{{ϵ}_1}}+{\left|\frac{2 y}{L_y}\right|}^{\frac{2}{{ϵ}_1}}\right)}^{\frac{{ϵ}_1}{{ϵ}_2}}+ {\left|\frac{2 z}{L_z}\right|}^{\frac{2}{{ϵ}_1}}=1 .$$

(3)

Equation (2) is a special case of (3) with the ratio of powers: ${ϵ}_1={ϵ}_2=2/g$. The moduli signs are omitted in equation (2), since, the shapes of LiFePO₄ crystallites can be described with only even g < 30 [14, 27] and their 0-th genus topology [28], i.e. excluding torus-type shapes.

Figure 1 shows anisotropic crystallites and columns $M_1$ along the 1st crystallographic direction. These forms are also used in the TOPAS Rietveld refinement [25]. Replacing x₁ in (2) with $\frac{M_1}{2}$ , we obtain for the length $M_1$

$$M_1=L_1{\left(1-{\left((2 x_2)/L_2\right)}^g -{\left((2 x_3)/L_3\right)}^g\right)}^{\frac{1}{g}}.$$

(4)

For ${\overline{D}}_1\left(M_1\right)$we get

$$D_{1 XRD}\equiv {\overline{D}}_{1 }\left(M_1\right) = \frac{1}{V}\int M_1 dV = \frac{1}{V}\iint {M_1}^2{dx}_2{dx}_3=K_n {\overline{L}}_{1 },$$

(5)

where V is the volume of the crystallite. The equality (5) is obtained for an ellipsoidal crystallite with orthogonal axes, and the coefficient $K_n$= 3/4 [7,9] results from normalizing to the average linear length of the crystallites. It will be seen below that in some cases the dependence of $D_{1 XRD}$ on only ${\overline{L}}_{1 }$is retained, while $K_n$ depends on the shape of the crystallites and the parameters of their size distribution. Variants with non-orthogonal axes are dealt with in detail in [9].

The double integration in (5) is carried out over the integration region Reg of the central section of the ellipsoid perpendicular to the 1st axis,

$$Reg= {\left({\left((2 x_2)/L_2\right)}^g+{\left((2 x_3)/L_3\right)}^g\right)}^{1/g}\le 1.$$

(6)

The advantage of using Lamé's superellipsoids is the ability to determine their volume and cross-sectional area perpendicular to column M₁ [26, 29], respectively, by using beta functions

$$V=\frac{1}{g^2}{L}_1{L}_2{L}_{3}B\left(\frac{1}{g},\frac{2}{g}+1\right)B\left(\frac{1}{g},\frac{1}{g}+1\right) ,$$

(7)

$$S_{\left(010\right)}=\frac{1}{g }{L}_2{L}_{3}B\left(\frac{1}{g},\frac{1}{g}+1\right)$$

(8)

In particular, knowing the mass of the cathode powder, its density and the statistics of the size distribution of crystallites, with the equations (7, 8) it is possible to determine the inner surface of the cathode and use it to improve the accuracy of determining the diffusion coefficient using the Randles-Sevcik equation [14], as well as for quantitative calculations of the rate of electrochemical reactions using the Butler–Volmer equation [27]. In this case, the shape of crystallites must be preliminarily determined by analyzing TEM and SEM images of crystallites either visually or using software processing methods [29-32].

2.2. Crystallite statistics

In the study of the structural and electrochemical properties of LiFePO₄ powders, Lognormal [33,34,35,36,37,38,39], Normal [40,41,42], and Weibull [36, 37] distributions were observed earlier. The latter was also used to describe carbon conglomerates of anode electrodes [39, 43]. Gamma and Lognormal distributions, as noted above [10, 11], were used in the analysis of specially tested CeO₂ powders. There are the following structural and technological justifications for their use:

• In [44], it was shown that during the crushing of rocks, a Lognormal distribution of products occurs, which is due to the multifactorial character of the crushing process itself;
• In [45], the Normal distribution arises as a result of aggregative growth and oriented attachment of nanocrystals;
• In [43], the Weibull distribution results from the restructuring of the electrode morphology during battery usage;
• Gamma is associated with exponential and Normal distributions. In [46], the Gamma distribution was observed at a constant rate of formation of nuclei near the supersaturation threshold of the solution and also during the formation of raindrops in [47,48].

In the SI section, these 4 functions are used to construct the one-dimensional histograms of size (diameter) distributions for 8000 virtual spherical particles. They can be used at the first stage of visual selection of the experimental distribution type.

It was noted above that only the Lognormal [13,14] and Normal [12] distributions can be used to describe anisotropic particles, taking into account their pairwise covariances between the sizes. For the remaining two, they can also be adjusted to Lognormal and Normal with smaller errors compared to those arising when covariances are not taken into account.

In [49], marginal size distributions of crystallites in LiFePO₄ powders were determined along all three crystallographic axes, and in [14] the parameters of the 3-dimensional lognormal function $f_{LM}\left(\overline{L}\right)$were determined

$$f_{LM}\left(\overline{L}\right)=\frac{1}{L_1{L}_2{L}_3\sqrt{{\left(2\pi \right)}^3\mathrm{d}\mathrm{e}\mathrm{t}\overline{\Lambda }}}\exp \left[-\frac{1}{2} {\left(\mathrm{l}\mathrm{n}\overline{L}-\mathrm{l}\mathrm{n}\overline{\overline{L}}\right)}^T{\overline{\Lambda }}^{-1}\left(\mathrm{l}\mathrm{n}\overline{L}-\mathrm{l}\mathrm{n}\overline{\overline{L}}\right)\right],$$

(9)

where $\overline{L}=\left[\begin{array}{c}{L}_1\\ L_2\\ L_3\end{array}\right]$ – crystallite sizes, $\overline{\overline{L}}=\left[\begin{array}{c}{\overline{L}}_1\\ {\overline{L}}_2\\ {\overline{L}}_3\end{array}\right]$ – their mean values, $\overline{\Lambda }=\left[\begin{array}{ccc}{W}_1^2& r_{12}{W}_1{W}_2& r_{13}{W}_1{W}_3\\ r_{21}{W}_2{W}_1& W_2^2& r_{23}{W}_2{W}_3\\ r_{31}{W}_3{W}_1& r_{32}{W}_3{W}_2& W_3^2\end{array}\right]$– matrix of correlation moments; the off-diagonal elements are covariances between the marginal distributions, for example, ${Cov}_{12}={Cov}_{21}=r_{12}{W}_1{W}_2$ between the 1st and 2nd, $r_{12}$ is the correlation coefficient between them. The possible values of $r_{12}$ range from 0 to 1.

Below we also use the Normal distribution $f_N\left(\overline{L}\right)$ in the form [15]

$$f_N\left(\overline{L}\right)=\frac{1}{\sqrt{{\left(2\pi \right)}^3\mathrm{d}\mathrm{e}\mathrm{t}\overline{\Lambda }}}\exp \left[-2 {\left(\overline{L}-\overline{\overline{L}}\right)}^T{\overline{\Lambda }}^{-1}\left(\overline{L}-\overline{\overline{L}}\right)\right].$$

(10)

The most significant difference between (10) and (9) is the linear scale used in (10), and the dimension of the dispersion: in (9) the dispersion W is dimensionless, and in (10) it has the dimension of length. In the calculations, it will be normalized to the average length along the i-th direction, which makes the calculations more general.

2.3. X-ray size calculation

According to (1), the X-ray diffraction size ${\overline{D}}_{1 }\left(M_1\right)$ along the 1st axis is the average of the lengths of the columns $M_1$, shown in Figure 1, over all powder crystallites. We obtain it by averaging over the i-th crystallite using integration (5) over the cross section (6) followed by averaging over all N powder crystallites using function (9) or (10). To do this, we use the Mathematica 12 software for element-by-element matrix multiplication and their summation in the following form, for example, for the Lognormal distribution

$${\overline{\overline{D}}}_1\left(M_1\right)=\mathbf{T}\mathbf{o}\mathbf{t}\mathbf{a}\mathbf{l}\left[\right(\overline{f_{LN}}\circ {\overline{D}}_{1 }\left(M_1\right) ],$$

(11)

where $L_{in}{, L}_{jn}$and$L_{kn}$are the sizes of crystallites along the [010], [100] and [001] axes, respectively, the index n varies from 1 to N, also for Normal $\overline{f_N}$. The crystallite volumes V_N in the equation (6) will be represented by a 3-dimensional matrix (7), depending on the g-type of superellipsoid, in particular for an ellipsoid with g=2, $V_N=\frac{\pi }{6} L_{in}{\ast L}_{jn}\ast L_{kn}$. $\overline{f}$ in (8) stands for the discretization of the function (9) or (10) in the form of a matrix normalized to 1, each element of which is the probability of superellipsoid presence with the corresponding size in the powder.

Thus, the task of numerical calculations (11) is to determine the dependence of the coefficient $K_n$in (5) on g – the type of superellipsoids and on the following parameters of the Lognormal (9) or Normal (10) distributions: ${\overline{L}}_{i }$ – average sizes, W_i – dispersion, and $r_{ik}$– correlation coefficients of the matrix $\overline{\Lambda }$.

3. Results of calculations and their discussion

In the presence of a large number of parameters, it is necessary to establish their hierarchy, evaluate their possible ranges, and also determine the number of bins (discreteness or dimension) of the histograms [50, 51], which is equal to the matrix dimensions used in (11) in our case. It has been verified that the bin number should be at least 15-20 for a small discretization error, comparable to the errors of powder structural studies (see SI).

3.1. Lognormal distribution

Parameter ranges. For high quality LiFePO₄ powders, the parameters of the function $f_{LN}\left(\overline{L}\right)$ fall in the following ranges: ${\overline{L}}_{1,2,3} ~$ 40÷500 nm, $W_{1,2,3}$~ 0.3÷0.6, r_12,13,23 ~ 0.4÷0.7 [27].

To explain the calculation algorithm, the test results for 8000 crystallites and a 20-bin distribution histogram by analogy with [27, 49] are shown in Figure 2. SI contains similar figures and numerical values of 6-bit matrices, which were used to check and control the stages of calculations and to determine the discretization error. From Figure 2b it can be seen that the histogram of X-ray diffraction size distribution shifts to the region of large sizes. This usually occurs with an increase in the dimensionality of the averages: from 1D to 3D. The linear dependence of $D_{1 XRD}$ on L₁ exactly corresponds to the value of K_n = ¾ (¾ of ${\overline{L}}_1$ ) in (5), Figure 2a (right scale). The sum of 270.41 nm is much larger than 240 nm, that suggests K_n > ¾, which is larger than expected for ellipsoidal crystallites with a small dispersion value W.

Dependences of K_n on correlation coefficients. The calculations illustrated in Figure 2 ignore the dependence of K_n on the correlation coefficients $r_{ij}$ and on their anisotropy. Figure 3 shows the results of testing these assumptions and draws the following conclusions:

• 1. After normalization ${\overline{\overline{D}}}_1\left({\overline{L}}_1\right)$ , the coefficients K_n do not depend on the values of the average sizes ${\overline{L}}_1$ in the tested range from 80 to 600 nm, i.e. the normalization significantly expands the universality of the developed calculation procedure.
• 2. A weak dependence of K_n on r is observed at r < 0.7, i.e. in the range of values observed in high quality LiFePO₄ samples. However, as shown in [27], this dependence is quite significant for the values of their electrochemical capacitances. Therefore, here and below the calculations are performed for the value r₁₂ = r₁₃ = r₂₃ = r = 0.5.
• 3. The observed dependence at r > 0.7 is related to the histogram discreteness. Only in this region anisotropy of crystallite sizes along other crystallographic directions can affect the calculation results. It also expands the universality of the developed calculation procedure.

Dependencies of K_n on the dispersion W₁. Figure 4 shows the results of calculations that can already be directly used in the development of technology for electrode crystalline powders with a higher accuracy in determining their structural parameters compared to those previously used in [14, 27] to determine the block structure of crystallites [49], to assess the potential of the developed technologies, and to optimize it [14]. From Figure 4, the following conclusions can be drawn:

• 1. For the Lognormal crystallite size distribution the following analytical expression for K_n is found

  K_n =( K_g,0 + K_g W² )

  (12)

  where K_g,0 is a constant at W→0, K_g is a constant depending on the g-type shape. The values of these constants are collected in Table 1.
• 2. Equation (12) is valid up to the values of W₁ ~ 0.6. For larger values, it is necessary to increase the bit length of the histograms, which is possible with a 2÷3-fold increase in the number of measured particles. On the other hand, the error of using 6-bit histograms with 2÷3 times smaller number of particles (up to 1000÷2000) will be no more than 1%.

3.2. Normal distribution

Parameter ranges. The examples of 20-bit histograms for 8000 particles given in SI show their main feature: at W > 1, the proportion of small-sized particles increases due to the impossibility of negative values of the particle size, and at W < 1, it is located entirely on the positive part of the abscissa axis. Thus, the average dimensions of ${\overline{L}}_{1,2,3}$ ~ 40–500 nm and the value of the normalized dispersion covering the point W=1 were used in the calculations.

Dependences of K_n on the correlation coefficients. Based on Figure 5, the following conclusions are made:

• At small W < 0.5-0.8, there is no dependence of K_n on r and W, as well as no dependence of the calculated values ${\overline{D}}_{1 }\left(M_1\right)$on the values of the distribution parameters along the 2nd 3rd axes, as is observed in Figure 3 for the values of r > 0.7 of Lognormal distributions. The values of K_n calculated at small W, up to nearly 1, can be considered universal constants equal to K_g,0 with a small error (Table 1).
• For large W up to 2, the dependences along the 2nd and 3rd axes are significantly weak. In particular, as can be seen from Figure 5 for the curves with W=2, the 30% changes in the values of dispersion lead to much smaller changes of less than 3% in K_n.

Thus, the calculations with equal parameters along all three axes can be considered accurate. With anisotropic deviations of up to 30% in the distribution parameters at large W, the error of 3% can occur.

Dependences of K_n on dispersion W₁. Figure 6 shows the results of calculations that can also be directly used in the development of technology for electrode crystalline powders with Normally distributed crystallite sizes to determine block structure [49], assess the potential of the developed technologies, and optimize it [14].

As can be seen from Figure 6, the solid lines for W₁=W₂ =W₃ can be approximated by two straight lines in the region W₁ < 2.4 that intersect at the point W = 1. Analytically, the situation can be described by the so-called piecewise-continuous function [52]

$$K_n=K_{g,0}\left\{\begin{array}{c}1 ;W<1 \\ \frac{2}{3}+\frac{1}{3}W ;W=1÷2\end{array}\right.,$$

(13

here the values of the constants $K_{g,0}$ are the same as in Table 1.

Obviously, the rather simple equations (12, 13) can be obtained as a result of analytical transformations within the framework of the corresponding approximations. However, a simpler way was to carry out numerical calculations with subsequent approximation of the calculated curves.

4. Conclusions

The use of the BWA (Bertaut-Warren-Averbach) method is an important tool for studying the structure of crystalline powders, namely, for relating the anisotropic X-ray diffraction sizes of crystallites to their anisotropic linear sizes (1). First, the lengths of the columns are averaged over the volume of a typical crystallite. Here the anisotropic shape of the crystallite is important. Second, the lengths of the columns are averaged over all crystallites of the powder under study. Here, the distribution of crystallites over their anisotropic linear dimensions is important. A sign of these averaging is the presence of a relation for the dimensions (length)⁴/(length)³. In particular, therefore, the result of averaging over a crystallite is divided by the volume of this crystallite. Using these principles of the BWA method, the following results were obtained:

• Numerical calculations have been carried out for the crystallites with orthorhombic lattice, such as LiFePO₄, NMC. The shape of crystallites was approximated with the Lamé's g-type superellipsoids including cuboids at large g > 20÷30. The values of the dimensional parameters of the powders were close to the ones observed experimentally.
• In the presence of anisotropy, the equations (12, 13) can be normalized to the length of the crystallite along its i-th axis; the normalized (dimensionless) distribution dispersions can also be used for Lognormal and Normal distributions. The universality of the calculation results is limited by the discreteness of the matrices, which is due to the discreteness of the experimental histograms of crystallite sizes in the powders under study.
• For Lognormal distributions, the coefficient K_n depend weakly on the correlation parameter r for r < 0.7, i.e. in the range of values observed in high quality LiFePO₄ samples. The observed dependence at r > 0.7 is related to the discreteness of the histograms and to the effect of crystallite size anisotropy along other crystallographic directions.
• For the Lognormal size distribution of crystallites, the following universal analytical expression for K_n (12) is found, which is applicable up to values of W₁ ~ 0.6. For larger values, it is necessary to increase the bit number of the histograms, which requires a 2-3-times increase in the number of measured particles.
• For Normal distributions with small W < 0.5-0.7, there is no dependence of K_n on r and W, and the calculated values ${\overline{D}}_{1 }\left(M_1\right)$are independent of the values of the distribution parameters along the 2nd 3rd axes. The calculated K_n for small W ≤ 1, can be considered as universal constants.
• For Normal distributions at large W ≤ 2, the dependences of K_n on the parameters of other axes are significantly weak, in particular, for the curves with W₁ = 2, 30% changes in dispersions along the 2nd and 3rd axes lead to much smaller changes in K_n of no more than 3%.
• For Normal distributions, the dependences of K_n on W for the ratios W₁=W₂ =W₃ can be approximated by two straight lines that intersect at the point of W = 1. Analytically, the situation can be described by the so-called piecewise-continuous function (13).
• For Gamma and Poisson distributions the mathematically reliable ways to describe correlations between sizes along different coordinate axes are lacking. In some cases, these distributions are visually close to the Lognormal and Normal functions which can be used instead of Gamma and Poisson.
• An elliptical cylinder (bar) (shown in Figure 1с) can be described as a combination of an ellipsoid with M₁ columns along the x₁ axis and a cuboid with M₂ columns along x₃.
• The use of the BWA method is carried out on the example of LiFePO₄ and is justified by the fact that this compound has been studied quite well, so the development stages and final conclusions can be varied. In the other case of crystalline powders with a non-orthorhombic lattice, it is necessary to take into account the angles between the axes [9]. In this case, it becomes possible to widely use the BWA method for studying various crystalline powders: other electrode and atmospheric contaminants [47], as well as nuclear fuel and materials [53,54].