The Effect of Crystal Defects on 3D High-Resolution Diffraction Peaks: A FFT-Based Method

Forward modeling of diffraction peaks is a potential way to compare the results of theoretical mechanical simulations and experimental X-ray diffraction (XRD) data recorded during in situ experiments. As the input data are the strain or displacement field within a representative volume of the material containing dislocations, a computer-aided efficient and accurate method to generate these fields is necessary. With this aim, a current and promising numerical method is based on the use of the fast Fourier transform (FFT)-based method. However, classic FFT-based methods present some numerical artifacts due to the Gibbs phenomenon or “aliasing” and to “voxelization” effects. Here, we propose several improvements: first, a consistent discrete Green operator to remove “aliasing” effects; and second, a method to minimize the voxelization artifacts generated by dislocation loops inclined with respect to the computational grid. Then, we show the effect of these improvements on theoretical diffraction peaks.


Introduction
X-ray diffraction (XRD) is one of the most powerful non-destructive tools to investigate materials, as their wavelength is commensurate with the distance between atoms within a crystal [1][2][3][4][5][6][7][8][9]. Successive improvements of both the X-ray sources (from X-ray tubes to third generation synchrotrons) and detectors (from photographic plates and gas counters to fast two-dimensional arrays) have led to a tremendous increase in the quantity of data recorded per unit time, allowing real time in situ or in operando measurements [10,11]. It is now possible to determine the 3D grain microstructure of a bulk material with a submicron resolution (using topo-tomography), to follow the evolution of the elastic strain state of the grains of a polycrystal during mechanical tests (3D-XRD, far field diffractometry), or to measure the distribution of strains within a few grains in real time (2D diffractometry) [12,13]. Such experiments result in terabytes of data recorded within a few days, which need to be analyzed efficiently. In fact, only a low fraction of those data is actually treated because scientists lack both time and numerical tools (software) for further analysis [14].
The classical techniques used to analyze the 1D or 2D diffraction patterns recorded during tests performed on polycrystalline specimens such as the Rietveld method, the square sines method to measure internal stresses, or CMWP (convolutional multiple whole profile) fitting for

FFT-Based Algorithm and Mechanical Fields
Let us consider a homogeneous elastic medium with eigenstrain assuming a periodic unit cell discretized in N × N × N voxels and subjected to a uniform overall strain tensor denoted E. Here, this overall strain is the spatial average of the strain field in the unit cell (with external loading and a given eigenstrain field). The unit cell may contain voxel size defects (0D) such as chemical inhomogeneities, line defects (1D) with arbitrary shape and distribution, planar (2D) defects such as stacking faults, or 3D precipitates which are all modeled with an eigenstrain tensor. (See [14].) These defects create a displacement field and thus generate strain and stress fields [38].
The displacement vector is denoted u and in the forthcoming equations, x denotes any position vector within the unit cell. All vector and tensor fields will be written using bold characters.
Starting from the equation for mechanical equilibrium, div σ(x) = 0, and using field equations (strain compatibility, generalized Hooke's law, decomposition of total strain compatible strain into elastic strain and eigenstrain), the displacement field is given at every position by the Green's function technique [39]: where the symbol * denotes the spatial convolution product, c 0 is the homogeneous linear elastic stiffness, ε * is the eigenstrain field and B is a third order Green operator defined in Fourier space as: in which B is the Fourier transform of B and G is the Fourier transform of the elastic Green tensor [39]. Therefore, using the Fourier transform of spatial convolution product, Equation (1) can be written in Fourier space as: Several numerical results showed that the use of the third order operator B derived from the classic Green G leads to spurious oscillations on the computed displacement field near materials discontinuities and dislocations [37]. The discrete Fourier transform (DFT) used in this algorithm indeed transforms a periodic function in real space into a periodic function in reciprocal space. However, the operator B commonly used is the continuous analytic operator truncated to the size of the unit cell of the reciprocal space: it is not periodic function. To fix this problem, we very recently developed a periodized consistent discrete Green operator using the DFT. The mathematical derivations of this discrete Green operator are given elsewhere [37]. The Fourier transform of this discrete Green operator denoted B is written as function of B and reads: Discrete frequencies appearing in this equation are given when N is even by (T is the period of the unit cell): Here, the sum on the B operator is extended to the whole reciprocal space (in practice for m, n, p up to a few tens) and folded up onto the unit cell of the DFT with suitable coefficients. The inverse transform of u(ξ) gives the displacement field at the center of each voxel.
We can also compute the displacement field at each voxel's corner with a shifted operator using the shift theorem:

Numerical Examples
Let us consider a homogeneous material with isotropic elastic constants: Young's modulus E = 333.4 GPa and Poisson ration ϑ = 0.26. This approximately corresponds to the room temperature elastic constants of single crystalline Ni-based superalloys. The unit cell (Figure 1a) is discretized in 128 × 128 × 128 voxels and contains a square-shaped inclusion discretized in 32 × 32 × 1 voxels corresponding to an Eshelby-like square prismatic loop perpendicular to the z-axis. In order to generate a shift of the upper surface of the inclusion relative to its lower surface by a Burgers vector b(0, 0, b 3 ), only the voxels within the inclusion are submitted to a non-zero eigenstrain tensor defined as: ε * ij = 0 except ε * 33 = 1. Then, we have b 3 = t × ε * 33 where t the thickness of the inclusion in the z-direction (i.e., the voxel size). This displacement field computed with the FFT algorithm using the different Green operators defined in Section 2.1 is represented along z-axis in Figure 1.

Voxelization Effect on the Displacement Field
While the displacement field computed for dislocation loops having their planes parallel to the faces of the simulated volume is correctly given with discrete Green operators ′ or ′′ , voxelization artifacts appear for inclined loops, as shown with ′ in Figure 2 for a dislocation loop with a [01 ̅ 1] Burgers vector lying in a (111) slip plane of a fcc crystal. The eigenstrain tensor is constrained in the region occupied by the dislocation loop (transformed voxels) and is given by: When computed along a line crossing a dislocation loop, the displacement field exhibits a discontinuity with a jump equal to Burgers vector b. This is indeed observed in Figure 1. However, the displacement field computed with the usual Green operator B (Figure 1b) also shows spurious oscillations as soon as the discontinuity is approached. These oscillations (numerical artifacts) are not observed with the periodized operators B and B . An artificial damping of the oscillations in Figure 1b (such as a low pass filtering) might smooth these oscillations, but it would also smooth the discontinuity, which is not searched.

Voxelization Effect on the Displacement Field
While the displacement field computed for dislocation loops having their planes parallel to the faces of the simulated volume is correctly given with discrete Green operators B or B , voxelization artifacts appear for inclined loops, as shown with B in Figure 2 for a dislocation loop with a 011 Burgers vector lying in a (111) slip plane of a fcc crystal. The eigenstrain tensor is constrained in the region occupied by the dislocation loop (transformed voxels) and is given by: where A s is the area on which planes with normal n(n 1, n 2 , n 3 ) has slipped by a relative amount b(b 1, b 2 , b 3 ) and V is the volume occupied by the dislocation loop [40,41]. As before, the dislocation loop is 32 voxels wide in the x and y directions, and 1 voxel thick but now with a z position such that x + y + z = constant. The displacement has been computed at the center of voxels with the periodized operator B along z ( Figure 2b). As in Figure 1c, the displacement in the center of a voxel belonging to the loop plane (black dot in the reddish transformed voxel in Figure 2b) is zero. The displacement in the first neighboring voxels (red dot in Figure 2b) are shifted relative to the expected position, so that the displacement difference between these voxels is significantly lower than b, see Figure 2c. It can be checked in Figure 2b that each of these voxels shares three faces with a transformed voxel. A more detailed analysis shows that the second neighbors (which share three edges with transformed voxels) are also slightly shifted in the opposite direction. The result is shown in Figure 2d with: a strong localized oscillation of the phase (taken here as the displacement modulo b).
Although the amplitude of this shift is small (less than 10% of the Burgers vector) it has unwanted consequences on the diffraction peak simulation: The dislocation loops are surrounded by four impaired layers of voxels: As the scattered X-ray amplitude is proportional to the Fourier transform of G·u (see Equation (8) in Section 4), we can expect a phantom streak in the intensity in a direction perpendicular to the loop plane.
The displacement field near the edges of the loop (near the dislocation line) will be quite different from its expected value, and the strain field will not vary with the distance r to the dislocation line as 1/r. This will strongly affect the tails of the diffraction peaks.
The dislocation loops are surrounded by four impaired layers of voxels: As the scattered X-ray amplitude is proportional to the Fourier transform of G· u (see Equation (8) in Section 4), we can expect a phantom streak in the intensity in a direction perpendicular to the loop plane.
The displacement field near the edges of the loop (near the dislocation line) will be quite different from its expected value, and the strain field will not vary with the distance to the dislocation line as 1⁄ . This will strongly affect the tails of the diffraction peaks.

Sub-Voxelization Method
The (conceptually) simplest way to remove this voxelization artifact would be to work on a multiple grid (to multiply the number of voxels along each direction by 2, 4, or more), then to downsample the displacement field data. In that case, FFT algorithms would lose much of their interest due to these more demanding computational efforts. We show below that this can be done in a more "economical"-and simple-way by applying a patch to the FFT-computed displacement field. The basic method is to compute, on the same grid, the difference vector: where ( ) is the displacement vector calculated for voxels where this eigenstrain is concentrated on a single plane of sub voxels ( Figure 3b) and ℎ the displacement field in direction i of voxels with a uniform eigenstrain (Figure 3a). For the sake of clarity, we use 2D diagrams in Figure 3, but here the technique is applied to real 3D problems.

Sub-Voxelization Method
The (conceptually) simplest way to remove this voxelization artifact would be to work on a multiple grid (to multiply the number of voxels along each direction by 2, 4, or more), then to downsample the displacement field data. In that case, FFT algorithms would lose much of their interest due to these more demanding computational efforts. We show below that this can be done in a more "economical"-and simple-way by applying a patch to the FFT-computed displacement field. The basic method is to compute, on the same grid, the difference vector: where u sub i (x) is the displacement vector calculated for voxels where this eigenstrain is concentrated on a single plane of sub voxels (Figure 3b) and u hom i the displacement field in direction i of voxels with a uniform eigenstrain (Figure 3a). For the sake of clarity, we use 2D diagrams in Figure 3, but here the technique is applied to real 3D problems.  In order to compute the displacement due to sub-voxels, we use a × ( × × ) grid for 2D (resp. 3D) problems where each voxel can be subdivided into × ( × × ) sub voxels. Only n ( × ) sub voxels are submitted to an eigenstrain field. At a point of the grid (black dots, Figure 4a), we need to compute the sum of the displacements due to the n ( × ) sub voxels j (center ) within a voxel centered at point .
This sum is equivalent to the sum of the displacements due to a strained sub voxel at point on the grid points such as = (Figure 4b). It is also equivalent to the sum of the displacements ′ due to a full voxel at point on the initial grid on points ′ such as ′ = (Figure 4c). The only difference between these last two sums is due to the long-range strain field, and approximately results in a linear drift of the displacement. As the end of the vectors ′ does not lie on the grid points (voxel centers) but on the corners of the voxels, the ′ displacements must be calculated with the shifted operator ′′ (Equation (5)). A last point is the scaling of the and ′ sums during the operations of Figure 4. To keep the one Burgers vector jump between both sides of the sub voxels plane in Figure 4a, the eigenstrain in the sub voxels must be multiplied by n. The backwards change of scale requires a division by n: there is no scaling factor between ℎ and = ∑ ′ .  In order to compute the displacement due to sub-voxels, we use a N × N (N × N × N) grid for 2D (resp. 3D) problems where each voxel can be subdivided into n × n (n × n × n) sub voxels. Only n (n × n) sub voxels are submitted to an eigenstrain field. At a point A of the grid (black dots, Figure 4a), we need to compute the sum of the displacements u j i due to the n (n × n) sub voxels j (center B j ) within a voxel centered at point O. This sum is equivalent to the sum of the displacements due to a strained sub voxel at point O on the grid points A j such as OA j = B j A (Figure 4b). It is also equivalent to the sum of the displacements u j i due to a full voxel at point O on the initial grid on points A j such as OA j = nB j A (Figure 4c). The only difference between these last two sums is due to the long-range strain field, and approximately results in a linear drift of the displacement. As the end of the vectors OA j does not lie on the grid points (voxel centers) but on the corners of the voxels, the u j i displacements must be calculated with the shifted operator B (Equation (5) two sums is due to the long-range strain field, and approximately results in a linear drift of the displacement. As the end of the vectors ′ does not lie on the grid points (voxel centers) but on the corners of the voxels, the ′ displacements must be calculated with the shifted operator ′′ (Equation (5)). A last point is the scaling of the and ′ sums during the operations of Figure 4. To keep the one Burgers vector jump between both sides of the sub voxels plane in Figure 4a, the eigenstrain in the sub voxels must be multiplied by n. The backwards change of scale requires a division by n: there is no scaling factor between ℎ and = ∑ ′ .  We need to compute ∆ ij_pl (x) the difference in displacement in direction i due to a voxel which belongs to the plane "pl" (for fcc "pl" is equal to (111), 111 , 111 , 111 ) of a dislocation loop with a Burgers vector j at a position x relative to the transformed voxel. In practice, in a material with cubic symmetry, it is sufficient to compute ∆ 13_(111) (x) and ∆ 33_(111) (x), and to use the symmetries of the cube (fourfold [001] axis, threefold [111] axis, and 110 symmetry plane) (and suitable exchanges of the components of x) to obtain the required components. As can be seen in Figure 2b, ∆ ij_pl (x) is non-zero only for the neighbors of the transformed voxel, except the drift due to the long range strain alluded above. The final computational procedure to determine ∆ ij_pl (x) and use the patch is now detailed: Compute the field c 0 : ε * defined in Equation (1) for an isolated voxel with the eigenstrain associated to a dislocation loop (Equation (6)) with a Burgers vector [001] in a (111) plane (see Figure 2).
Compute the displacement field in directions x (u hom 1 ) and z (u hom

3
) at the voxels' center around the transformed voxel by convolution with the discrete periodized operator B (Equation (4)).
Compute the displacement field in directions x and z at the voxels' corners around the transformed voxel by convolution with the shifted operator B (Equation (5)).
Use the farthest voxels to correct the drift of the components so that all terms for large x are zero, and keep non zero only the terms for the first three neighbors.
The patch can then be applied on the raw (FFT-based) displacement field by adding the convolution of all transformed voxels of the different slip systems by the relevant ∆ ij_pl (x).

Results
For numerical tests, we used the same 128 × 128 × 128 grid as before, and the transformed voxel was divided into 8 × 8 × 8 sub voxels (using a reference medium with the same elastic constants as before). Only the final values in units of b (after drift correction) of ∆ 13_(111) (x) and ∆ 33_(111) (x) are used and other components are obtained by symmetries (Appendix A).
The patch was used on the same configuration as in Figure 2. Figure 5a shows the resulting displacement field and Figure 5b the phase (i.e., the displacement modulo a Burgers vector) in Burgers vector units. As it can be observed from Figure 5a, the voxelization artifacts of the displacement field are removed thanks to the sub-voxelization method. In addition, the resulting phase varies smoothly even during the crossing of the dislocation loop which is more realistic, see Figure 5b.
Use the farthest voxels to correct the drift of the components so that all terms for large are zero, and keep non zero only the terms for the first three neighbors.
The patch can then be applied on the raw (FFT-based) displacement field by adding the convolution of all transformed voxels of the different slip systems by the relevant ∆ _ ( ).

Results
For numerical tests, we used the same 128 × 128 × 128 grid as before, and the transformed voxel was divided into 8 × 8 × 8 sub voxels (using a reference medium with the same elastic constants as before). Only the final values in units of b (after drift correction) of ∆ 13_(111) ( ) and ∆ 33_(111) ( ) are used and other components are obtained by symmetries (Appendix A).
The patch was used on the same configuration as in Figure 2. Figure 5a shows the resulting displacement field and Figure 5b the phase (i.e., the displacement modulo a Burgers vector) in Burgers vector units. As it can be observed from Figure 5a, the voxelization artifacts of the displacement field are removed thanks to the subvoxelization method. In addition, the resulting phase varies smoothly even during the crossing of the dislocation loop which is more realistic, see Figure 5b.

Application on Diffraction Peak Simulation
In this section, we show simulated diffraction peaks in order to point the effects of voxelization artifacts and of the patch on numerical results. Under kinematical conditions and assuming a coherent beam, the amplitude of a diffracted wave at a position q in the vicinity of a reciprocal G lattice vector is [14,18,[42][43][44]: where x is the position of the scattering atom, A 0 (x) is the amplitude of the incidence wave, F(G, x) is the local structure factor, and u(x) the displacement field. The scattered intensity is I(q) = |A(q)| 2 . For a face-centered cubic crystal, this intensity is non zero when G (h, k, l) is such as h, k, and l have the same parity. Here two diffraction vectors G (200) and G (002) are used. They respectively correspond to G·b = 0 and G·b = 1. The 3D diffracted intensity has been calculated using the FFT instead of the continuous Fourier transform, then summed in the planes perpendicular to the G vector to obtain a linear plot along G equivalent to a I(2θ) plot. In Figure 6a (G (200)) and Figure 6c (G(002)), we show the diffracted intensity (logarithmic scale) as a function of the pixel position i, and in Figure 6b,d a logarithmic/logarithmic plot of the intensity vs. |i − i 0 | where i 0 is the center of the peak. In order to only study the effect of the displacement fields, we set A 0 (x) = 1 and F(G, x) = 1 for these simulations.
diffracted intensity has been calculated using the FFT instead of the continuous Fourier transform, then summed in the planes perpendicular to the G vector to obtain a linear plot along G equivalent to a I(2θ) plot. In Figure 6a (G (200)) and 6c (G (002)), we show the diffracted intensity (logarithmic scale) as a function of the pixel position , and in figures 6b and 6d a logarithmic/logarithmic plot of the intensity vs. | − 0 | where 0 is the center of the peak. In order to only study the effect of the displacement fields, we set A 0 ( ) = 1 and ( , ) = 1 for these simulations.  The peak shape near the top of the peaks is the same for both computing methods. It is perfectly symmetric in the G·b = 0 case and exhibits a bump on the right side for G·b = 1. The long-range behavior is, however, quite different. When the displacement field has been calculated with the usual truncated operator (black line), a phantom peak is observed at large |i − i 0 | (at large q), which is due to the short period oscillations near the displacement field discontinuity (Figure 1a). The behavior of the peak calculated with the modified Green operator (red curve) is only slightly better: the intensity at large q is underestimated in one case and overestimated in the other. When the intensity has been calculated with the sub voxel patch (dark blue curve) the long range intensity follows the expected I 0 |i − i 0 | −3 law [45,46]: the peak tails are indeed related to the highly distorted zones near the dislocations' cores. However, the dark blue curve saturates at very large q. We suppose this is due to the use of the FFT instead of the continuous Fourier transform in the calculation of the scattered amplitude (Equation (8)). The plot of Figure 6 represents only one period in Fourier space, and is repeated over and over on all Fourier space. We can now calculate the intensity of the tails of these repetitions: where m varies from −5 to 5 (zero excluded). If we now plot the difference between the dark blue intensity curve and this background line, we obtain the pink curve. On the log./log. plots, Figure 6b,d, it can be checked that this curve follows the I 0 |i − i 0 | −3 law to the end. Thus, the residual error in the intensity computed by FFT results of the FFT itself, and not from an error on the sub voxel-corrected displacement field. If the number of voxels is increased to 512 3 or 1024 3 while keeping the physical size of the representative volume constant, this residual error should fall down to undetectable levels.

Conclusions
In this paper, we have shown that although the use of a periodized Green operator in the FFT-based method improves the final displacement field solution in a representative volume containing discontinuities (dislocation loops), artifacts due to the voxelization of the dislocation loop planes are still present with respect to analytical solutions. These artifacts have unwanted consequences on the tails of diffraction peaks simulated by using this displacement field as input data.
We have introduced a patch through a sub-voxelization method, which corrects these artifacts by simulating the displacement field without employing a finer grid resolution. A simple construction method for this patch has been given and the patch can be used in a single post-processing step to modify the initial FFT-based displacement field.
The modified displacement field has been used to simulate one-dimensional diffraction peaks. The procedure strongly improves the shape of the peaks' tails, i.e., it gives a good description of the displacement field and the phase near the dislocation lines.