Vibrational Properties of Nanocrystals from the Debye Scattering Equation

One hundred years after the original formulation by Petrus J.W. Debije (aka Peter Debye), the Debye Scattering Equation (DSE) is still the most accurate expression to model the diffraction pattern from nanoparticle systems. A major limitation in the original form of the DSE is that it refers to a static domain, so that including thermal disorder usually requires rescaling the equation by a Debye-Waller thermal factor. The last is taken from the traditional diffraction theory developed in Reciprocal Space (RS), which is opposed to the atomistic paradigm of the DSE, usually referred to as Direct Space (DS) approach. Besides being a hybrid of DS and RS expressions, rescaling the DSE by the Debye-Waller factor is an approximation which completely misses the contribution of Temperature Diffuse Scattering (TDS). The present work proposes a solution to include thermal effects coherently with the atomistic approach of the DSE. A deeper insight into the vibrational dynamics of nanostructured materials can be obtained with few changes with respect to the standard formulation of the DSE, providing information on the correlated displacement of vibrating atoms.


Atomistic simulations
The Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS, [5]) code was employed to simulate the dynamics of 16,727 Palladium atoms arranged in a fcc (lattice parameter a = 3.89Å) sphere of 19 unit cells (7.40nm) diameter. Atomic interactions were ruled by the Embedded Atom Method (EAM, [6,7], with interatomic potential from [8]) and the system was evolved in the microcanonical (NVE) ensemble using 1fs timestep at 300.00 ± 0.13K. The same procedure was reiterated to the most representative aggregate of the set employed to model scattering data, the particle depicted in figure M4b and 1b, composed of 197,587 atoms (the cube edge was 38 unit cells or 15.057nm and T = 300.07 ± 1.01K).

Time-average scattering data
The time-averaged pattern was computed by averaging the outputs of the DSE (equation M1) applied to 500 snapshots of the trajectory separated by 1ps. The (spherical) coherent scattering factor was Correlation of thermal displacements can be characterized both evaluating the projection of the displacement of i-th atom on that of atom j, δ i (t) ·δ j (t), and the pair thermal displacement, with r ij (t) = r i (t) − r j (t) and r ij = ⟨r ij ⟩ = ⟨r i ⟩ − ⟨r j ⟩. The dimensionless coefficient is also defined, with ⟨u 2 ⟩ being the average Mean Squared Displacement (MSD). It must be noticed that being the above averages computed in a three-dimensional space whereas the quantity in equation M8 resides in a two-dimensional space, a conversion factor is needed to link atomistic and scattering quantities, e.g. for the MSD The average number of atomic pairs separated by a distance r at time t is given by the Pair Distribution Function (PDF aka radial pair correlation function, see e.g. [10]), being δ (x) the Dirac-delta function [11], N the number of atoms composing the aggregate and ρ the numeral density. The PDF was applied to the same set of frames used to compute the time-averaged scattering pattern therefore obtaining a time-averaged pair distribution function. Coherently with the enforced harmonic approximation, to provide a different view and compare results of motion correlation from scattering data analysis, PDF peaks were fitted using Gaussian functions (see e.g. [12] and figure 2). Given the variance σ 2 ij , the relation allows to compare the different views.

Scattering data modeling
Both simulated and experimental scattering data were modeled using the Stochastic Real-space Modeler (StoRM, [13]) code, implementing a simulated annealing algorithm [14][15][16]. The quantity to be minimized is the difference between intensity observed (I o ) and computed for a given model (I m ), being σ the standard uncertainty, P the number of pattern points and ν the varied parameters. Following, the most important aspects when dealing with experimental data are reported. Starting from a careful analysis of High-Resolution Transmission Electron Microscopy (HRTEM) pictures, sixteen (atomistic) particles were built according to crystallographic principles. The form factor associated to each point scatterer was the same described by equation 1.
Background modeling Air scattering and Kapton ® contribution were accounted for by measuring an empty capillary. Data were fitted with a mixture of mathematical functions, namely seven pseudo-Voigt curves and a fourth-order Chebyshev polynomial (definitions are in [13]). This fit function was then perturbed with a third-order Chebyshev polynomial of the first kind during the modeling of data emanated from Palladium nanocubes and added to the output of the TDSE. The two curves are reported in figure 3, to highlight the background contribution.

Deformation of crystallographic particles
The position of the i-th atom in the deformed configuration can be expressed as [13,17,18] The fit parameter a is responsible for a constant offset of the deformation which insists on the particle, normalized for convenience to the bulk lattice parameter a 0 , whereas σ modulates a term due to Pauling [19,20], not null for undercoordinated atoms, i.e. atoms with less than N = 12 (for the fcc system) nearest neighbors. The function κ produces a hkl-dependent behavior, being sr the compliance of the crystal projected on directionr, s its average value over all possible crystal directions and α a fit parameter. Considering cubic crystal symmetry, the reciprocal of the Young's modulus alongr = {u 1 , u 2 , u 3 } can be expressed in terms of compliance matrix (S) elements as (see e.g. [21]), The average atomistic strain can then be evaluated by computing the difference of bond length (the Euclidean norm of r ij = r i − r j , i being a nearest neighbor of j) and averaging the quantity over the number of nearest neighbors, The average strain of atomic bonds of a given atom n as seen by its nearest neighbors can therefore be expressed as Varied parameters To summarize, 5 parameters were fitted to account for the background (1 for scaling the fit function and 4 for a third-order Chebyshev polynomial), 16 to define each particle fraction and 5 for the deformation model. Indeed, both values of a and σ were mapped on different sizes through a Young-Laplace-like law [22,23], being e the edge of the cube and ϕ either a or σ. The number of parameters for vibrational properties ranged from 2 (⟨δ 2 ⟩ and Θ for the Debye correlated model) to the chosen number of k ij values to be refined plus one (⟨δ 2 ⟩).
Sensitivity to k ij parameters For a given set of varied parameters, data modeling was reiterated ten times with different seeds for the (ranlux [24]) pseudo-random number generator, implemented within the GNU Scientific Library [25]. To assess the effect of the number of fittable k ij parameters on residuals, i.e. the difference between intensity observed and computed for a given model, equation 8 was evaluated for each configuration. Figure 4a illustrates the evolution of χ 2 when progressively increasing the number of k ij parameters different from unity together with the fraction of neighbors in a given shell, which reveals the sensitivity to each parameter, as demonstrated from a different point of view in figure M5.
To conclude the analysis, figure 4b illustrates the output of equation M9 applied to the particle depicted in figure 1, both implementing the uncorrelated model (k ij = 1 for each shell) and using fifteen k ij parameters from experimental data modeling. parameters, from zero (corresponding to the uncorrelated model since the relation k ij = 1 holds for all parameters) to fifteen. Yellow bars depict the population of a given coordination shell (the most representative particle has been analyzed) and justify the relative weight of the pattern component associated to a given k ij (see also figure M5b). Right, simulation of the scattering pattern emanated from the particle sketched in figure  1 assuming uncorrelated motion (k ij = 1 for each shell) and with fifteen k ij parameters from experimental data modeling.

Statistical analysis
For the complete set of varied parameter, a typical χ 2 curve is reported as a function of the annealing temperature in figure 5a. Collecting only combinations of parameters giving a χ 2 reasonably close to the "best value", in this case values in [ min{χ 2 }, 1.05 min{χ 2 } ] , i.e. inside the red region depicted in figure 5a, an unbiased variance-covariance matrix can be drawn by computing each coefficient as The diagonal of the variance-covariance matrix defines the variance of a given entry n, reported in figure 5b for a selected set of parameters as coefficient of variation The Pearson product-moment correlation coefficients [26], for a typical case are also reported in figure 6, limiting the analysis to the parameters describing the background (ι and a 0 to a 3 ), the particle fraction (ϕ 1 to ϕ 16 ) and atomic vibrations (⟨δ 2 ⟩ and k (1) ij to k (15) ij ). Interestingly, while background and particle fraction terms exhibit strong correlation, vibrational coefficients (especially ⟨δ 2 ⟩ and k ij values for low index neighbor shells) are rather uncorrelated from every parameter.

Longitudinal and transversal vibration modes Equation M9
can be generalized to allow for independent transversal and longitudinal (with respect to the bond axis) vibrations by relieving the condition ⟨sin 2 θ ′ ⟩ = ⟨cos 2 θ ′ ⟩ = 1/2. Therefore, k ij parameters for longitudinal (k (L) ij ) and transversal modes (k (T ) ij ) can be introduced in equation M7, leading to where,

Series expansion of equation M4 with terms up to O
can be easily made more accurate by adding higher order terms of the series expansion of equation M4, leading to being cos (Qr ij ) } .
The above expressions, including terms up to O ( ⟨δ 2 ⟩/2 ) 3 , are sufficiently accurate for the proposed case of study. However, adding further terms to the expansion of equation M4 is relatively straightforward, although it involves an increasing algebraic complexity.

Additional material
To highlight the effect of atomic vibrations on the powder diffraction pattern emanated from the most representative particle of the set employed to model experimental data, animation 1 illustrates the effect of k ij values associated to the first three neighbor shells. Additionally, the line corresponding to uncorrelated motion (all the k ij parameters equal to unity, green line) and the pattern corresponding to a static object (⟨δ 2 ⟩ = 0Å 2 , blue line), implying the TDSE to reduce to the DSE, are also reported.
Animation 1: Correlation of atomic vibrations: Effect of k (1) ij to k (3) ij parameters on the powder diffraction pattern for the most representative particle of the set employed to model experimental data. If atoms are fixed, ⟨δ 2 ⟩ = 0Å 2 and the TDSE (equation M9) reduces to the DSE (blue line). The line profile associated to atoms vibrating independently (uncorrelated motion) and according to a Gaussian distribution is depicted by the green line (⟨δ 2 ⟩ = 0.0156Å 2 and every k ij equals unity). Given the same ⟨δ 2 ⟩ of the previous case, the red line simulates different combinations of k (1) ij , k (2) ij and k