Accurate quantification of lattice temperature dynamics from ultrafast electron diffraction of single-crystal films using dynamical scattering simulations

In ultrafast electron diffraction (UED) experiments, accurate retrieval of time-resolved structural parameters, such as atomic coordinates and thermal displacement parameters, requires an accurate scattering model. Unfortunately, kinematical models are often inaccurate even for relativistic electron probes, especially for dense, oriented single crystals where strong channeling and multiple scattering effects are present. This article introduces and demonstrates dynamical scattering models tailored for quantitative analysis of UED experiments performed on single-crystal films. As a case study, we examine ultrafast laser heating of single-crystal gold films. Comparison of kinematical and dynamical models reveals the strong effects of dynamical scattering within nm-scale films and their dependence on sample topography and probe kinetic energy. Applying to UED experiments on an 11 nm thick film using 750 keV electron probe pulses, the dynamical models provide a tenfold improvement over a comparable kinematical model in matching the measured UED patterns. Also, the retrieved lattice temperature rise is in very good agreement with predictions based on previously measured optical constants of gold, whereas fitting the Debye–Waller factor retrieves values that are more than three times lower. Altogether, these results show the importance of a dynamical scattering theory for quantitative analysis of UED and demonstrate models that can be practically applied to single-crystal materials and heterostructures.


I. INTRODUCTION
Ultrafast electron diffraction (UED) has emerged as a powerful tool for structural dynamics research, allowing one to record excitedstate atomic structure evolution with sub-picosecond temporal resolution. [1][2][3][4][5] Recently, UED has been applied to study diverse phenomena, including nonequilibrium phases and transformations in quantum materials, 6 formation of warm dense matter, 7 and electron-phonon coupling mechanisms, 8 to name a few. The broad utility of the technique lies in the sensitivity of diffraction signals to many structural features such as crystal structure, lattice strain, atomic coordinates, and thermal displacement parameters. 9 While analysis of UED data has often examined the evolution of diffraction signals themselves, additional important information can be gained by quantitatively retrieving the underlying structural parameters. For instance, time-resolved lattice temperature informs the role of heat flow in the observed dynamics and can be used to determine other quantities such as electron-lattice coupling constants, 10,11 time constants for defect and interfacial scattering, 12 and thermal conductivity within or between layers. 13 Also, accurate time-dependent structure retrieval may reveal new transient phases or finer structural details of metastable phases.
However, many solid-state UED studies are performed on dense, oriented single-crystal films, 1,[6][7][8]14 which pose challenges for quantitative structural retrieval. Single crystals are often chosen, because they provide numerous advantages for UED, including strong peaks associated with individual diffraction orders, access to diffuse scattering between the peaks to study phonon population dynamics, 8,15,16 lack of grain boundary scattering which factors into polycrystalline film dynamics, 17 and opportunity for polarization-dependent study and control. 18 In addition, some emerging materials are primarily available as single crystals. 19,20 On the other hand, the high scattering cross section for electron probes often leads to multiple scattering, and in single crystals, this is compounded by electron channeling down atomic columns. 21 These effects are especially strong when probing dense, inorganic solids along high-symmetry zone axes. Despite the reduced cross section for relativistic electron probes, such effects have still been observed in some experiments even at MeV-scale energies. 14 When these effects dominate, kinematical approximations are no longer valid, and modeling the complete "dynamical" diffraction process is required to match the diffraction signals.
Suitable dynamical scattering models for UED of single-crystal films are needed for accurate quantification, but efforts to develop such models have so far been limited. A multiple scattering theory has been applied to UED of crystal surfaces in the reflection geometry 22 and to strain wave imaging in UEM. 23 As for UED in the transmission geometry, a few studies have invoked Bloch wave eigenvalue (also called "N-beam") [24][25][26][27] or multislice 28 simulations to improve fitting to measured signals as well as accuracy of retrieved lattice temperature and phonon dynamics. However, many of these examples added empirical functions to quantitatively match the diffraction signals, indicating room for further improvement in the underlying models. For instance, some factors that are not always considered for TEM simulations become important for UED, including lattice temperature, partial coherence of the probe, and sample topography (averaged over the typical lm to mm UED probe size).
In this article, we demonstrate dynamical scattering models that are suitable for matching UED signals from single-crystal films and retrieving the lattice temperature dynamics. We first describe the computational approaches used, including both a multislice and a Bloch wave method, and introduce adaptations to account for key physical parameters. We then illustrate the role of dynamical scattering in UED of single-crystal films by comparing static and temperature-dependent diffraction signals calculated using kinematical and dynamical models for gold films of varying thicknesses and rippling as well as varying electron probe energy. Finally, we apply these models to analyze relativistic UED measurements of single-crystal gold films recorded at the High Repetition-rate Electron Scattering (HiRES) beamline of Lawrence Berkeley National Laboratory. 29 We show quantitative matching of static UED patterns, obtaining a factor of ten improvement from the dynamical models over the kinematical model and achieving an R factor of 2%. We then demonstrate lattice temperature retrieval, showing that dynamical scattering models provide good agreement with expectations based on the known optical properties of gold, whereas the kinematical model underestimates the expected temperature rise by nearly three times.

II. UED SIMULATION METHODS
The simulation approach in this work was developed to account for several important factors in UED experiments of single-crystal foils, which are illustrated in Fig. 1(a). An excitation, such as a laser FIG. 1. Modeling ultrafast electron diffraction signals from single crystal foils using the dynamical scattering theory. (a) Illustration of a UED experiment on a freestanding crystal foil, highlighting some key features to include for accurate models. (b) Channeling plots showing an example of strong dynamical scattering effects. These were simulated using the mutlislice method for a 750 keV electron wave through gold [001] with mean square displacements of u 2 ¼ 0:024 Å 2 . The amplitude, jwj (relative to incident amplitude of 1), and phase, / (in radians), across a single gold unit cell through a 40 nm thick crystal are shown. (c) Graphical summary of the procedure for UED simulations performed in this work. Images in step (i) are of a projected potential (V 2D p ) slice for gold [001] before and after applying Debye-Waller (DW) damping to account for thermal displacements: u 2 ¼ 0:024 Å 2 (T ¼ 300 K) is shown as an example. Note these images are each normalized by the maximum value for clearer visualization. Image in step (ii) shows a generated 2D Gaussian orientation distribution, which is sampled using increasingly dense grids like shown by the pink dots until convergence. Diagram in step (iii) illustrates computation of a thickness-dependent stack of diffraction patterns at each orientation. Diffraction patterns for beam-sample angles of 0 mrad (left), 50 mrad (center), and 100 mrad (right) at 1 (top), 6 (middle), and 11 (bottom) unit cells are shown as examples.
pulse, triggers a dynamic process, leading to numerous potential excited states. The diffraction signals for each state are an average over the lm to mm-scale probed region, which can consist of a wide beam-sample orientation distribution due to the sample topography and the probe beam divergence. Also, the electron probe can experience multiple scattering and channeling effects as described by the dynamical scattering theory: for example, Fig. 1(b) shows a simulation of the amplitude and phase of the envelope of a 750 keV electron wave propagating through a 40 nm thick gold crystal oriented along [001]. Within just a few nm, the amplitude and phase become highly non-uniform and show complex variation with thickness.
Our procedure is summarized in Fig. 1(c). For each proposed ground and excited state, we generate the electrostatic potentials for the crystal, generate possible beam-sample orientation distributions, simulate the thickness-dependent diffraction patterns for the sampled orientations, and then compute weighted sums of the patterns according to the proposed distribution. In this section, we will first discuss the underlying scattering models used and then describe how the beam-sample orientation distributions and excited states (in this case, crystals with varying lattice temperatures) were incorporated. We note here that the Bloch wave and multislice methods give nearly equivalent results over the parameter ranges studied, so we use these methods interchangeably throughout the article (see the Appendix).
A. Models for diffraction from a single-crystal film

Kinematical scattering
Formulas for calculating diffraction signals in the kinematical approximation have been described in detail elsewhere, 9,30 so we only elaborate on the details specific to our approach here. In this work, the "weak phase object" (also called "Moliere") approximation is used to compute the atomic scattering factors, f e;j , from parameterized atomic electrostatic potentials, V a , as calculated by Kirkland using a relativistic Hartree-Fock program 30,31 f e ðqÞ ¼ 2pi k where J 0 is the zeroth order Bessel function, q is the reciprocal space distance, and r is the real space distance. We use the relativistic electron interaction parameter, r e , defined as for de Broglie wavelength k, kinetic energy E 0 , electron rest mass m 0 , speed of light c, and electron charge e. We note that we did not use absorptive electrostatic potentials in this work because, in our case, they mostly remove electrons that should in fact remain included. The dominant contribution to absorptive potentials is typically thermal diffuse scattering (TDS). 32 However, the angle spread of the UED probe in our case is large enough that much of the TDS remains in the measured diffraction peaks and should not be removed from the simulations. Future work could examine modifications to absorptive potentials for such cases, though in our case we expect it would only provide minor corrections.
The scattering factors are then used to compute the structure factor for the periodic crystal, F hkl , for each diffracted beam at reciprocal lattice vector g hkl with Miller indices h, k, and l Applying the shape factor for a thin film gives the diffracted intensity as a function of film thickness, t, 33 where s hkl is the excitation error and n hkl is the extinction distance where V cell is the unit cell volume and b is the angle between the beam and the surface normal. As shown in Fig. 1(b), the weak phase object approximation that underlies this method can be rapidly violated in dense, oriented crystals even at relativistic beam energies. Along atomic columns, the phase is strongly disturbed, leading also to strong modification of the amplitude envelope. This will be examined further in Sec. III.

Bloch waves
For thicker specimens where kinematical approximations no longer hold, scattering patterns can instead be calculated by solving the Schr€ odinger equation for the electron wave passing through the specimen. The Bloch wave eigenvalue solution is convenient for crystals. The electron wave and specimen potential are decomposed into Fourier components, and a matrix equation is derived by which the electron wave components (and, hence, diffracted intensities) can be computed at varying distances through the crystal.
The derivation of the matrix equation and the approximations used here are given in Ref. 30. The computational procedure is to first calculate the Fourier components of the scattering potential, U hkl , where h is Planck's constant, F Ã hkl is the complex conjugate of F hkl calculated using the first Born approximation, 30 and the other symbols are defined above. In addition, the excitation errors s hkl are calculated.
Then, a subset of the Fourier components is selected to be included in the simulation: namely, those with non-zero U hkl , in-plane reciprocal space distance q xy < q xy;max , and s hkl < s max . The thresholds q xy;max and s max are set such that the diffracted beam intensities of interest are converged (see Subsection II D). For this work, q xy;max ¼ 4:5 Å À1 and s max ¼ 0:1 Å À1 . Note that Fourier components beyond those of the signals of interest are included since the diffracted beams interact with each other.
The scattering potential components and excitation errors are used to build the matrix A in which where k 0 is the incident electron wave vector and i, j are indices for the Fourier components included in the simulation. Notably, computing Structural Dynamics ARTICLE scitation.org/journal/sdy U g j Àg i typically requires computing U at scattering vectors outside of the Fourier components selected for the simulation. The eigenvalues, 2k 0;z c j , and eigenvectors, C j , of this A matrix can then be used to compute the electron wave w in terms of the chosen Fourier components as a function of depth, z, where C is a matrix with eigenvectors C j as the columns. The diffracted beam intensities are then

Multislice
Another solution to the aforementioned Schr€ odinger equation is dividing the specimen into a series of 2D projected potential slices, interacting the wave with each slice, and then propagating to the next slice. This multislice approach is widely used for electron microscopy image simulation and has been described in detail elsewhere, 30,34 so it will only be described briefly here.
In this work, the envelope of the electron wave function, w, is initialized as a plane wave. It is then advanced through each slice j with thickness Dt by applying two operators sequentially. First, the interaction operator is applied in real space Here, V 2D p ðrÞ is the projected electrostatic potential within the slice, computed using the parameterized atomic potentials as determined by Kirkland,30 and r e is the relativistic interaction parameter as defined in Eq. (2). Second, a propagation operator is applied in reciprocal space Finally, the diffraction signals are obtained from the Fourier transform of the exiting envelope Since the simulations in this work are oriented along the [001] zone axis of a cubic crystal, we choose the slices to be equally spaced, each containing one layer of atoms. Also, the simulation cell is a single unit cell in the models used here. We note that modeling thermal diffuse scattering, such as by using the frozen phonon approach, requires larger simulation cells. 35 The image size was 256 Â 256 px, chosen to achieve convergence (see Sec. II D).

B. Incorporating thermal motions
For this work, the lattice temperature is incorporated by applying Debye-Waller damping to the projected potentials. This approximation models the electron beam traveling through a time-averaged electrostatic potential and has been shown to account for the influence of thermal motions on the coherent Bragg diffraction peaks. 35 As explained in Sec. II A, we expect thermal diffuse scattering (TDS) to have only minor effects on the measured peak intensities in our case and so it is not included here. Where TDS plays an important role, diffraction simulations with finer q resolution over multiple "frozen phonon" configurations could be averaged together, 35 though at much greater computational expense.
Atoms moving randomly and independently with an RMS displacement of u RMS along each dimension form a 2D Gaussian distribution of positions in the plane This distribution is convolved with the electrostatic potential in real space, effectively acting as a Gaussian filter. This filter in Fourier space is This approach is readily generalized to anisotropic or anharmonic thermal motions by applying the corresponding two-dimensional filter.
In the kinematical theory, this filter is applied to the diffracted intensities as a Debye-Waller factor (DWF) When kinematical approximations are valid, this allows extraction of a change in RMS displacements, Du RMS , from measurements of diffraction intensities for a set of diffraction orders through linear regression using the form Meanwhile, in the Bloch wave approach, this filter is applied to the scattering potential components On the contrary, in the multislice approach, this filter is applied to the projected potential slices An example of this is shown for the first slice from a gold [001] unit cell at a temperature of 300 K in Fig. 1(c-i).
In the dynamical scattering models, the relationship between Du RMS and changes in diffracted intensities is more complex and depends greatly on both intrinsic and extrinsic sample properties such as material, orientation, thickness, and sample topography, as will be illustrated in Sec. III. A simple analytical formula like Eq. (15) does not generally exist for dynamical scattering; instead, Du RMS can be determined by minimizing the least squares error between simulated and measured diffraction intensity changes.

C. Orientation averaging
In many UED experiments, a distribution of beam-sample orientations is sampled simultaneously due to two factors: sample rippling Structural Dynamics ARTICLE scitation.org/journal/sdy within the large (often mm-scale) lateral probe size and angular spread of the beam due to partial coherence. We model this by computing and incoherently summing the diffraction signals over a distribution of tilt angles, i.e., by computing the following weighted integral over the 2D orientation space, A, where h x and h y are the horizontal and vertical tilt angles, respectively, and pðh x ; h y Þ is the distribution of orientations. A similar approach has been used to model precession electron diffraction (PED), 36 and this approach has been employed by others to account for UED beam divergence. 37 Notably, this is a 2D integral and has significant contributions from near-zone orientations, so it will be more expensive to compute and more sensitive to the sampling than 1D integrals for PED which may avoid sampling near the zone axis. In these calculations, pðh x ; h y Þ is approximated as a circularly symmetric Gaussian distribution, like illustrated in Fig. 1(c-ii). This distribution can be considered as the convolution of the beam angular spread with the sample rippling such that the RMS tilt spread, r h , is given by r h;probe is typically on the order of 1 mrad or less for transmission UED setups, 37 and r h;sample depends on the sample preparation but can be as much as tens of mrad in some cases. 14 Since the calculations in this work are for a zone axis with fourfold symmetry, gold [001], the sampling grid for a given r h just spans the positive quadrant in orientation space from 0 to 3r h , and then fourfold rotational averaging is applied to account for the other quadrants. Points located more than 3r h from the center are set to zero to maintain circular symmetry. In this way, %99% of the Gaussian volume is sampled.
We compute the integral over the orientation space using an iterative 2D trapezoidal quadrature algorithm. 38 On the first iteration, a square sampling grid is initialized with just four samples: one at each corner. Then, on each successive iteration, points are added to complete a sampling grid with half the spacing in each dimension. A running integral is computed by adding 3/4 of the newly integrated points to 1/4 of the previous integral value. The sampling points for iterations 3 and 4 of this procedure are shown in Fig. 1(c-ii) as an example.
The computations in this work sampled a 480 Â 480 mrad 2 tilt range to examine r h up to 160 mrad, though the plots only show up to 120 mrad for clearer visualization. For most simulations presented here, nine iterations were used for the entire tilt range sampled (giving 1.86 mrad sample spacing), and an additional iteration was performed for the inner 120 Â 120 mrad 2 (giving 0.94 mrad sample spacing). The temperature-dependent library used for Sec. IV B was only computed up to 20 nm film thickness, so a 3.72 mrad sample spacing for the whole range and 1.86 mrad spacing for the inner quarter was sufficient.
This algorithm is robust to the complex variations in the diffraction intensities with tilt angle. Especially convenient is the hierarchical nature of this algorithm: Each sampling grid can also be used to calculate tilt-averaged diffraction patterns for smaller tilt spreads, i.e., the grid used to compute the iteration N for r h is the same to compute the iteration N -1 for r h 2 . This property allows us to compute a library of tilt-averaged patterns with varying r h largely in parallel with additional iterations applied to the successively smaller tilt spreads once the larger tilt spread calculations are converged. For kinematical and Bloch wave simulations, the sample tilt is incorporated into the excitation error coefficients s hkl . For the Bloch wave simulations, this also affects which Fourier components are included in the simulation at each tilt angle, as different beams are brought near their Bragg condition and contribute to the scattering process.
For the multislice calculations, we implement sample tilt by applying the Fourier shear theorem to the propagation operator, as follows: 30 This allows one to sample arbitrary tilt angles without changing the electrostatic potential slices. Prior works have suggested that this approximation can introduce significant error at angles beyond 1 , 30 while freestanding films studied in UED sometimes have more than 5 RMS tilt spread. So, while this approximation is invoked here, future works could seek to implement approaches that improve accuracy for larger tilt spreads. 39,40

D. Convergence
To achieve good quantitative precision, some simulation parameters needed to be tuned until the diffraction signals converge. We used the crystallographic R factor as the metric, given by which compares iteration i and the final iteration, i max . All simulations in this work are converged until the R factor computed over the first seven diffracted orders (the ones we are interested in quantifying) is less than 1%. For Bloch wave calculations, the parameters to converge are the thresholds that determine which diffracted beams are included, i.e., the q xy;max and s max used here. For multislice calculations, the main parameter to converge is the real-space pixel size (q range), which must be sufficiently small (large) to include enough of the scattered beams for accurate diffraction calculations. The image dimensions are chosen to be powers of 2 for optimal speed of the fast Fourier transforms. Additionally, for all orientation-averaged simulations, the sampling density in the orientation space was increased until R < 1% was achieved over the entire thickness and RMS tilt range studied.

E. Programs
We have created a MATLAB code library to perform all the calculations shown in this work, which is available online (see Data Availability). All simulation methods can be performed using CPUs, but for the multislice calculations in this work, the time cost was prohibitive when performing the orientation-averaged calculations. As such, for the multislice calculations, we have also implemented a version for GPU computing using "gpuArray" objects in MATLAB.

ARTICLE
scitation.org/journal/sdy Our main goal with these codes was to demonstrate and compare the accuracy of the described approaches for the UED quantification shown, so speed was not fully optimized. For reference, each orientation-averaged diffraction library (thickness from 0 to 40 nm and r h from 0 to 160 mrad) took about 2.5 h to compute using our multislice program on an NVidia Quadro K5000 GPU and about 1.7 h using our Bloch wave program on an Intel iCore i7-8550U CPU. This time was reduced to about 23 min for the libraries computed up to just 20 nm thickness for Sec. IV B using Bloch waves. For all methods, the initial setup computations required several seconds, and then each orientation took less than a second to compute. Several approaches to improve performance have been demonstrated by others and could be implemented in the future. Examples for Bloch waves have included using off diagonal matrix elements to compute an array of tilts simultaneously 41 and GPU acceleration. 42 We also note that open-source, high-performance programs for Bloch wave and multislice simulations have been developed by others, 34,43 though some adaptation may be needed to suit UED simulation.

III. ROLE OF DYNAMICAL SCATTERING IN UED OF SINGLE-CRYSTAL FOILS A. Static diffraction peak signals in flat and rippled foils
The importance of dynamical scattering in oriented single-crystal foils is evident in the complex evolution of the electron wave through nanometer-scale ultrathin films even at relativistic beam energies. Here, we show diffraction simulations computed for a 750 keV electron wave passing through gold films up to 40 nm thick oriented along [001] with mean square displacements of u 2 ¼ 0:024 Å 2 . (This u 2 has been measured by others for films at 300 K using x-ray diffraction.) 44,45 Figure 1(b) shows channeling plots calculated using the multislice method, which illustrate the amplitude and phase of an electron plane wave passing through a flat film. Whereas kinematical scattering approximations assume the material is a weak phase object, the strong and dense gold atomic columns locally shift the electron phase by more than p within just a few unit cells. This imparts dramatic modifications in the electron wave amplitude, including significant channeling along the atomic columns within a few nm. This in turn leads to complex, oscillatory thickness dependence of the primary and diffracted beam intensities as shown in Fig. 2(a). Within just a few nm, the intensities deviate from those predicted by kinematical theory (dashed lines) despite the relativistic electron beam energy.
When averaging over a large beam-sample orientation distribution, much of the oscillatory behavior is smoothed out; however, the diffracted intensities still show a complex behavior as a function of thickness that is not captured by the kinematical theory. Figure 2(b) shows the total diffracted intensities calculated for a Gaussian tilt distribution with a large r h ¼ 100 mrad as might be found in a strongly rippled thin-film sample, which nonetheless shows significant variation in the diffracted intensities within this thickness range.
The deviation of dynamical scattering models from the kinematical theory can be quantified by computing the crystallographic R factor between the diffracted intensities obtained using two methods (I hkl;1 and I hkl;2 ) where the scaling factor a is fit to minimize R. The R factors computed between kinematical and dynamical simulations, R KinÀDyn , computed using the first seven diffracted orders for films of varying thickness and tilt spread are shown in Fig. 2(c). Larger tilt spread increases sampling away from the zone axis, where channeling and multiple scattering effects are reduced, and averaging over a broad range smooths out these effects. Thus, larger tilt spread increases the range of film thickness for which the diffracted intensities can be approximately calculated using the kinematical theory. Still, in this case, R KinÀDyn > 10% is observed for films thicker than 8 nm even at large, %100 mrad RMS tilt spreads. In many UED experiments, RMS tilt spreads can be much smaller, especially if films are prepared on sturdy membrane supports, and deviations from the kinematical theory are more pronounced. 14 Notably, the large variations observed occur well within the elastic mean free path, calculated using the classic formula like in Ref. 46 Structural Dynamics ARTICLE scitation.org/journal/sdy to be 18.8 nm for 750 keV electrons through gold. However, the mean free path essentially considers an average electrostatic potential where scattering events are uncorrelated, while an oriented single crystal presents highly correlated scattering events along the atomic columns which more rapidly lead to strong multiple scattering and violations of kinematical approximations. This observation highlights the need to apply dynamical scattering models to UED signals from single-crystal foils at thicknesses well below the elastic mean free path.

B. Dependence on electron probe energy
One of the motivations for developing and utilizing UED beamlines with higher electron probe energy is to reduce dynamical scattering effects like those shown in Sec. III A. Here, we examine how electron probe energy affects the validity of kinematical approximations for the case of single-crystal gold films. To do this, we performed multislice simulations for the same film thicknesses and beam-sample orientation distributions chosen previously, but now for various probe energies between 30 keV and 4 MeV. We then extracted the film thickness at which R KinÀDyn first exceeds 10% for each r h and beam energy. We note that the 10% value was chosen as an example, and the true acceptable range for the kinematical approximation depends on many factors, including the structural parameters being quantified and the desired accuracy. Still, the results for selected r h shown in Fig. 3 highlight some interesting trends and provide a sense of scale.
Notably, we observe different behaviors depending on the beamsample orientation distribution. For flat foils and low-divergence probes (r h near zero), the acceptable thickness range gradually increases with the beam energy, though still limited to less than 3 nm even at 4 MeV. For highly rippled foils where r h is tens to hundreds of mrad, we observe an initial increase in the acceptable thickness range and then a plateau at about 7.5 nm beyond a threshold energy. The threshold energy appears to decrease for broader orientation distributions. This behavior can be rationalized as follows. In the flat case, strong oscillations dominate [see Fig. 2(a)] with periods set by the extinction distances, which increase with beam energy even at relativistic energies [analogous to the ones defined in Eq. (5)]. Meanwhile, orientation averaging smooths out this oscillatory behavior [see Fig. 2(b)], and the deviation from kinematical validity is instead set by the interaction parameter r e , which levels off for relativistic beam energies. 30 Altogether, these results show that using higher probe energies and orientation averaging can help to extend the validity of kinematical approximations, but only up to a point; in many cases, dynamical scattering models will still be required to accurately match the signals.

C. Temperature dependence of diffraction signals
Dynamical scattering not only affects the individual peak intensities but also how they change with structural parameters such as the lattice temperature. Figure 4 compares the diffraction peak intensity changes computed for a 105 K temperature increase (Du 2 ¼ 0.008 Å 2 ) using Bloch wave and kinematical scattering calculations. Whereas the kinematical theory predicts the peak intensities will be scaled by the Debye-Waller factor regardless of crystal thickness, Bloch wave calculations show significant deviations from this prediction within few-nm film thicknesses. This is especially apparent for a flat crystal oriented on zone [ Fig. 4(a)] where diffracted intensities can be nearly eliminated or dramatically enhanced (up to and exceeding a factor of 10) with the temperature change depending on the film thickness. Orientation averaging again smooths the variations, but there are still significant deviations within several nm thicknesses.
Importantly, as the film thickness in the rippled case increases, the mean absolute intensity change due to temperature tends to zero with some diffraction peaks gaining intensity and others losing intensity as the temperature increases in films thicker than %10 nm. This behavior is in stark contrast to the kinematical theory, where all diffraction peaks lose intensity with increasing temperature. This can be understood by considering that in the regime of strong multiple scattering, electrons are scattering back and forth between various diffracted beams and the primary beam, so an increase in temperature merely modifies the distribution of these multiply scattered electrons throughout the various beams. Indeed in the flat film, this even leads to conditions where diffracted beams are observed to increase in intensity on average with increasing temperature.
As a result, fitting the Debye-Waller factor (DWF) to quantify lattice temperature changes in the dynamical scattering regime can give large errors, even if fitting several diffraction orders. The error in DT extracted by a least squares DWF fit [Eq. (15)] to the intensity changes for the first seven diffraction orders calculated using the Bloch wave model for DT ¼ 105 K is shown in Fig. 4(c). Again, orientation averaging reduces oscillatory behaviors and smooths out the error, extending the range of validity for DWF fitting compared to flat films, but significant errors still emerge within several nm. Notably, DWF fitting tends to significantly underestimate DT in these rippled film models due to the trend toward zero mean intensity change observed in panel (b), even extracting nearly zero temperature change in rippled films near 20 nm thickness (error % À100%). The deviations worsen for smaller tilt spreads with DWF fitting massively over-or underestimating the temperature rise, at some thicknesses even extracting a temperature decrease instead of an increase (error < À100%).
The fitting results can also vary dramatically depending on which diffraction orders are included in the analysis and how they are weighted.  3. Dependence of the onset of dynamical scattering effects in single-crystal gold films on electron probe energy. The film thickness where R KinÀDyn , the R factor computed between kinematical and multislice calculations of the first seven diffraction orders, first exceeds 10% is plotted as a metric for four RMS tilt spreads (r h ).

ARTICLE
scitation.org/journal/sdy Altogether, these simulations illustrate the important role of dynamical scattering in UED of flat and rippled single-crystal foils even at relativistic beam energies using gold as an example. They also reinforce that kinematical scattering models are insufficient for quantitative matching and analysis in films with strong multiple scattering. In Sec. IV, we demonstrate that the dynamical scattering models shown here can be used in practice to quantitatively match experimental UED data and retrieve structural parameters.

IV. QUANTIFICATION OF PHOTOINDUCED LATTICE TEMPERATURE RISE IN SINGLE-CRYSTAL GOLD FILMS
Here, we apply the described scattering models to quantitatively match and analyze a UED experiment performed on a single-crystal gold foil at the HiRES beamline. A 750 keV electron probe was used with a 150 lm RMS spot size. The [001]-oriented freestanding single-crystal gold foil on a 3 mm diameter, 300 mesh TEM grid was purchased from Ted Pella and was quoted to be 11 nm thick.

A. Matching an experimental UED pattern
We first applied the dynamical scattering models to match a UED pattern recorded from the sample without any optical excitation. An optical micrograph of the film is shown in Fig. 5(a). Large rippling is evident in the freestanding gold foil. The experimental UED pattern recorded at HiRES is shown in Fig. 5(b). The peak positions in reciprocal space are consistent with those expected for the [001] orientation of gold. The apparent fourfold symmetry of the pattern indicates the film is well oriented along the zone axis on average and suggests that the orientation distribution can be reasonably approximated as isotropic. We also note the peak widths are dominated by the electron probe divergence r h;probe , measured to be about 0.33 mrad. In this case, the total r h of the beam-sample orientation distribution will be dominated by the large sample rippling.
All experimentally recorded diffraction patterns examined in this work were obtained by the following method. 30 sub-frames were recorded for each pattern, and a dark current reference (recorded with both the photocathode and pump lasers off) was subtracted from them. Then, the alpha-trimmed mean of the sub-frames was computed, removing outlier x-ray spikes. The peak intensities for the first  (c) R factor between measured and simulated signals (R ExpÀSim ) using Bloch waves (dynamical) and kinematical models over a range of crystal thickness and RMS tilt spread (r h ) with u 2 ¼ 0.024 Å 2 . The green circle marks the best fit using the dynamical model, where R ExpÀSim ¼ 2:1%. At the same thickness and tilt spread, the kinematical theory gives R ExpÀSim ¼ 21%.

ARTICLE
scitation.org/journal/sdy seven diffraction orders were extracted by fitting radially symmetric 2D Gaussians of the form where q hkl is the peak location in reciprocal space, I hkl is the peak intensity, k hkl is the (uniform) background level beneath the peak, and r q is the peak width. I hkl is the parameter of interest to fit for each peak in each pattern. q hkl is also refined for all peaks in each pattern to accommodate drift and thermal expansion, and k hkl is likewise refined to accommodate background fluctuations and underlying thermal diffuse background. Meanwhile, r q is fixed to the 200 peaks in the first laser off pattern, then fixed to this value for all peaks in all subsequent patterns; this is done because the peak width in our experiment is dominated by the angular spread of the probe. Only those peaks with a corresponding Friedel pair visible in the pattern were included to reduce error from slight misorientation relative to the zone axis. Agreement between the experimental intensities and the diffraction patterns simulated with u 2 ¼ 0:024 Å 2 for varying film thicknesses and tilt spread was quantified by computing R ExpÀSim [Eq. (22)] using both dynamical and kinematical models, displayed in Fig. 5(c). Remarkably, the dynamical scattering models achieve a tenfold reduction in R ExpÀSim , yielding 2.1% at the optimal thickness and RMS tilt spread compared to 21% obtained with kinematical models. Furthermore, the optimum parameters are physically reasonable: a thickness of 13.5 nm is in good agreement with the 11 nm quoted by the vendor, and the large 95 mrad RMS tilt spread is reasonable given the optically visible wrinkling. Both models are only weakly dependent on tilt spread beyond % 20 mrad RMS, so in this range, the precise value of tilt spread is less important; on the other hand, the dynamical scattering model provides a precise determination of the film thickness that the kinematical model cannot.

B. Matching photoinduced difference patterns
Next, we applied the simulations to retrieve the photoinduced lattice temperature from a pump-probe UED measurement. Photoexcited UED patterns from the same single-crystal gold film were measured using k ¼ 1030 nm pump laser pulses at a 0.5 kHz repetition rate with varying fluence. At each fluence, UED patterns were recorded as the pump-probe delay, Dt, was scanned from À17.3 to þ56.0 ps using 6.67 ps steps. A coarse sampling was used in these measurements with a focus on extracting fluence-dependent temperature rise rather than the fine temporal dynamics. The average difference pattern recorded after the arrival of laser pulses with 6.3 mJ cm À2 (from þ22.5 to þ56 ps) is shown in Fig. 6(a) as an example. The coherent Bragg diffraction peaks are generally suppressed, and the diffuse scattering background generally increases as expected for an increase in incoherent thermal motions. However, the diffraction peak intensity changes deviate from the scaling of the Debye-Waller factor: For instance, the 200, 400, and 600 peaks show little change while the 220, 420, and 620 peaks show large changes.
The time-dependent diffraction peak intensities were extracted from the UED datasets with Gaussian peak fitting, and average peak intensities before and after Dt ¼ 0 were calculated. The change in mean square atomic displacements was then retrieved before and after time zero both by fitting the Debye-Waller factor [Eq. (15)] and the Bloch wave models to the intensity changes measured for the first seven diffracted orders. In the Bloch wave approach, the thickness and RMS tilt spread of the film were fixed, and Du 2 was optimized by interpolating the peak intensities between calculations performed at 15 values of u 2 ranging from 0.024 to 0.038 Å 2 . Both kinematical and Bloch wave models are fit by minimizing the mean square error of Àlog I hkl I 0;hkl for the first seven diffracted orders.
The results of this procedure for the after time zero diffraction intensities for a peak fluence of 6.3 mJ cm À2 are illustrated in Fig. 6(b). Indeed, the intensity changes predicted using Bloch waves (green diamonds) are a better match to the observations (black circles) than are those predicted by the Debye-Waller factor (dashed line), reducing the least squares residual by about a factor of 3. The variations in the intensity change between orders are largely captured by the dynamical scattering models used here, though differences still remain, perhaps due to differences between the simulated and actual orientation distribution of the sample or to finer details not yet accounted for such as inelastic scattering effects. Photoinduced strain could also alter the orientation distribution of the sample and contribute to peak intensity changes; that said, we note that including Dr h as a second fit parameter (done in a separate analysis) does not have a strong systematic FIG. 6. Quantifying light-induced lattice heating from pump-probe UED measurements. (a) Example of a photoinduced difference pattern recorded at HiRES using a peak laser fluence of 6.3 mJ cm À2 . (b) Measured ratio of laser-on to laser-off diffraction intensities I on =I off (black dots), Debye-Waller factor fit (orange line), and Bloch wave model fit (green diamonds). (c) Extracted fluence-dependent changes in total RMS atomic displacements (D u 2 ) and lattice temperature (DT) using Debye-Waller factor fitting (kinematical) and Bloch wave calculations (dynamical). Changes predicted using known optical constants of gold are superimposed as a black dashed line for comparison.

ARTICLE
scitation.org/journal/sdy effect on the results with Dr h less than 3 mrad and Du 2 modified by about 610% on average for all the peak fluences studied. The photoinduced change in mean square displacements and lattice temperature rise are plotted as a function of the peak fluence in Fig. 6(c). The relationship between mean square displacements and lattice temperature rise in gold is approximately linear in the studied range of lattice temperatures (about 3.947 Â 10 4 K/Å 2 ). 44,45 Strikingly, the dynamical scattering models retrieve photoinduced lattice temperatures that are more than three times higher than those retrieved using the Debye-Waller factor approach.
Comparing to estimations of lattice temperature rise for the given peak fluence using the known optical constants of gold supports the accuracy of the dynamical scattering models. The lattice temperature rise DT ¼ T f À T i was calculated by relating the absorbed energy density (i.e., in J/mol) U abs to the heat capacity of the material C p , where F inc is the incident laser fluence, A is the absorbance of the material at the incident photon energy, V mol is the molar volume, and t is the thickness of the film. The absorbance in the 13.5 nm film of k ¼ 1030 nm light, using n ¼ 0.153 and k ¼ 6.654, 47 was calculated using the coherent transfer matrix method 48 to be 3.7%. Using this and the measured temperature-dependent heat capacity of gold, 49 the temperature rise per unit of incident laser fluence was found to be 11.6 K/(mJ cm À2 ). Fitting a line to the temperatures retrieved with Bloch waves gives a slope close to this of 12.3 K/(mJ cm À2 ) whereas the kinematical approach gives 4.0 K/(mJ cm À2 ).

V. CONCLUSIONS
This work demonstrated the importance and application of the dynamical scattering theory for quantitative analysis of ultrafast electron diffraction patterns. As shown here for single-crystal gold foils, diffraction signals are influenced by film thickness, temperature, and topography in ways that are sometimes entirely opposite of intuitions from the kinematical theory. By virtue of the proposed modified treatment, we are able to reach accurate UED pattern matching and lattice temperature quantification in a single-crystal experiment. We also show how a kinematical approach to the same problem would lead to greatly underestimated lattice temperatures.
The described models can be further extended to a wide range of experiments and samples. For instance, they are readily extended to multilayered single-crystal films by simulating each layer in series and can be extended for crystals of any space group. Larger simulation cells and complex symmetries may demand higher performance programs and computing resources for practical computation and refinements. Other physical parameters not included here may also be important for accurate quantification for different classes of specimen and different experimental setups. Some possible examples include anisotropic thermal displacement parameters, thermal diffuse scattering, core losses, and coherence properties of the probe.
Other structural parameters besides the lattice temperature can be refined. For instance, crystallographic parameters like Wyckoff positions can be used as variables, and refining these in combination with the thermal displacements could be used to quantify and separate simultaneous crystal structure change and thermal motions during structural phase transformations.
More widespread availability and use of dynamical scattering models for UED will enable more detailed information to be retrieved from UED experiments and will expand the technique's capabilities and scientific breadth. A few examples of materials that could become more accessible include thick, multilayered crystals such as large epitaxial superlattices; single-crystal nanowires or nanoparticles with critical dimensions in the dynamical scattering regime; and buried layers in thick semiconductor device stacks. In the long term, improving quantitative matching of UED patterns could ultimately enable full crystal structure refinement of transient structures such as photoinduced nonequilibrium phases. To recover complete 3D movies of the atomic coordinates and thermal motions in single crystals, for instance, from UED tilt series, would mark a major milestone for UED and provide detailed structural knowledge of transient intermediates, metastable phases, coherent lattice responses, and the overall energy flow and structural dynamics.

DATA AVAILABILITY
The data that support the findings of this study are openly available in GitHub at https://github.com/dbdurham/QuantUEDSim, Ref. 50.

APPENDIX: COMPARISON OF BLOCH WAVE AND MULTISLICE METHODS
In principle, Bloch wave and multislice approaches can be equivalent since they both are methods of solving the same Schr€ odinger equation for fast electrons through the specimen under nearly the same approximations. In practice, however, which method is more convenient, faster to compute, or more accurate depends on the specimen material and geometry, experimental conditions, and the signals being modeled. In our case of modeling flat and rippled single-crystal gold Bragg diffraction peaks, we find that the approaches give nearly identical results over most of the thickness and tilt spread range studied. The R factor between the diffraction signals calculated using Bloch wave and multislice models, R BWÀMS , is shown in Fig. 7. We find R BWÀMS < 1% for most of the range studied, though errors slightly increase for r h > 100 mrad perhaps due to limitations of the approximate tilt correction to the multislice calculation discussed in Sec. II C.
For the small unit cell and low q resolution of these simulations, we find our Bloch wave model is several times faster than our multislice model. For more complex, larger simulation cells with higher q resolution, multislice may become the faster approach due to more favorable scaling with the number of reciprocal space points. 30