GARFIELD, a toolkit for interpreting ultrafast electron diffraction data of imperfect quasi-single crystals Open Access

Crystal structure determination by single-crystal x-ray or electron diffraction is usually based on the collection of a series of diffraction images, each containing a number of sharp and distinct reflections. Typically, the first step in diffraction data analysis is some form of “interpretation” of the diffraction patterns, which leads to an enormous reduction in the amount of data. The goal is to identify isolated reflection peaks and assign them to reciprocal lattice points (RLPs) of a perfect crystal with Laue indices hkl and further to use the measured (partial) intensities to obtain estimates of the underlying full Bragg intensities. These full intensities are—in the simplest cases—proportional to $|Fhkl|2$⁠, the squared moduli of the structure factors. The indexing task can be demanding, but powerful auto-indexing tools have been developed for various use cases and incorporated into commonly used crystallographic software packages.^1–16 They all amount to some form of spatial or geometric analysis of tabulated peak positions. For large unit cells, a single diffraction pattern can be sufficient to determine the orientation of the diffracting crystal and uniquely assign the reflection indices, especially with prior knowledge of the unit cell. This is of great benefit for serial femtosecond crystallography (SFX) and serial synchrotron crystallography (SSX) of macromolecular crystals where thousands of diffraction images, each taken from a different crystal, can be processed in a short time and with a high success rate.^17–20 

Due to their small wavelength and strong interaction with matter, the diffraction of electrons opens up new possibilities, but it also presents new challenges. The contribution of multiple scattering to the diffraction signal can be significant, potentially necessitating a dynamical approach rather than a kinematical one. The samples must be sufficiently thin, typically on the order of 100 nanometers or less, depending on the atomic composition. Moreover, the intrinsic small wavelength renders the indexing of individual diffraction images more difficult, as each image only encompasses information from a nearly planar section through reciprocal space. To address these challenges, a range of innovative electron diffraction protocols have been developed in modern electron microscopy. The most crucial aspect is the limitation of the probe volume to thin and flawless regions of the sample through the use of selected area electron diffraction (SAED) or nanobeam diffraction (NBD). In recent years, notable advancements have been made in the field of crystal structure determination by “three-dimensional electron diffraction,”²¹ both in material science and macromolecular crystallography. The advent of new methods and automated protocols, including microcrystal electron diffraction (MicroED),²² precession electron diffraction (PED),²³ and parallel NBD combined with STEM imaging, has considerably alleviated the adverse influences of radiation damage as well as multiple scattering and facilitated indexing. The continuous expansion of the range of single-crystal diffraction data that can be indexed automatically has made serial electron crystallography with nanocrystalline materials²⁴ and protein nanocrystals²⁵ a feasible undertaking.

Ultrafast electron diffraction of bright femtosecond electron pulses (UED)^26–34 for the analysis of structural dynamics poses unique requirements for the interpretation of the diffraction data. The goal of a UED experiment is to study structural changes after perturbation of the sample, e.g., by a short laser pulse at time $t0$⁠, by recording a series of diffraction images at times $ti=t0+Δti$ and comparing them with images recorded without perturbation (⁠ $Δti<0$⁠). The data acquisition must be at a sufficient signal-to-noise level to allow accurate measurement of differences in diffraction intensities due to potentially small structural changes. UED data collected from thin, single-crystal samples typically consist of one or several time series of diffraction images, each series comprising a large number of nearly identical diffraction patterns. Thus, “indexing” one representative image per time series would be sufficient to index all members of the series. The challenges of interpreting such UED data often begin with the problem of assigning the correct Laue indices to the observed reflections since the observed data often do not look sufficiently similar to expectations derived from oversimplified approaches that are ultimately unsuitable for describing this type of diffraction data.

The generation of temporally short and laterally confined high-fluence electron pulses inevitably leads to space charge effects that ultimately limit beam coherence.³⁵ In order to mitigate these effects, the electron beam size at the sample in UED is usually quite large, typically on the order of 100 μm in diameter. The combined requirement of a small sample thickness for high electron transmission and a large area to overlap with the beam dimensions makes UED samples susceptible to bending, cracking, corrugation, and other effects that degrade the sample quality and thus the quality of the diffraction data. Hence, most UED-ready crystalline samples are far from perfect single crystals and are better described as an ensemble of crystalline domains or crystallites with similar but not identical orientation. In addition, the preparation of UED samples, e.g., by thin-slicing crystals, often results in a high density of lattice defects, which reduces the average size of coherently diffracting domains and increases the effective size of reciprocal lattice points (RLPs) in all three spatial directions. For these technical and physical reasons, the Bragg reflections recorded in UED experiments are typically blurred and often overlap with other reflections, so that purely geometry-based indexing, which relies on precise peak positions, cannot be regularly applied. Typical UED data are of lower quality if judged from the viewpoint of well-established techniques of crystal structure analysis by x-ray or electron diffraction. An example from UED is shown in Fig. 1(a). For the determination of static crystal structures, such images are usually excluded from further consideration at an early stage of the standard workflows, e.g., in the “crystal screening” phase of single-crystal structure analysis or by fast and efficient filtering and restriction of indexing to “hits” in SFX and SSX. However, despite their shortcomings, time series of UED images can provide valuable information on structural dynamics, making the interpretation of such data highly desirable. Fortunately, the number of diffraction images to be indexed is usually small, namely, one per time series. Therefore, a different approach is required to overcome the challenges of indexing typical UED data.

The Garfield software package, which is the subject of this report, has been developed with the objective of facilitating the step of indexing single-crystal UED data under challenging conditions. In essence, the software simulates electron diffraction within a model that is computationally efficient and fast while retaining the relevant features found in typical UED diffraction patterns. The software provides tools for determining model parameters (including orientation of the crystal lattice) that, according to reasonable criteria, best describe the observed diffraction data. It differs from available indexing tools in several aspects. The most distinct differences are as follows: (1) the use of reflection intensities and peak positions in the analysis of diffraction data and (2) the use of a graphical user interface allowing direct and interactive control of the evaluation process by breaking out of the common data-in-result-out black-box approach. The Garfield approach to interpreting and indexing diffraction data is an interactive process. The objective is to identify the “best” solution, which is approached step by step. This involves excluding unlikely solutions and optimizing promising solutions through non-linear least squares (NLS) fitting of model parameters that describe diffraction geometry and sample properties. The two core tools in this process are GridScan, performing fast grid searches of possible orientations in three-dimensional rotation space, and GeoFit, a flexible tool for NLS fitting of individual parameters and whole parameter sets with presently up to 22 parameters for a single diffraction pattern. These include the (mean) orientation of the crystal lattice, the beam center, intensity and spatial scaling, sample thickness, and the orientation distribution of crystal domains.

The paper is structured as follows: Sec. II formulates the two guiding principles used to design the Garfield software package and describes how these principles are reflected in the chosen approach for indexing typical single-crystal UED data. Section III summarizes the technical requirements for installing Garfield (A), provides some general instructions and recommendations for working with the graphical user interface (B, C), and introduces the core tools, GridScan (D) and GeoFit (E). Section IV describes the underlying model used to predict the position and intensity of diffraction spots for a given set of model parameters. In Sec. V, the essential characteristics of Garfield are summarized, its benefits and limitations are discussed, and ways to improve its performance are outlined.

Although Garfield was developed with UED applications in mind and tested with data from low-symmetry molecular crystals, it can be useful in a broader range of applications where conventional methods for indexing fail. Additional examples of the application of Garfield can be found in the supplementary material.⁷³ 

Indexing single-crystal diffraction data is primarily a task of unraveling the geometric relationship between the probe beam direction, crystal lattice orientation, and detector geometry, with challenges posed by imperfect beam and crystal conditions coupled with detector insufficiencies. Conventional approaches rely on the precise determination of the positions of the recorded Bragg peaks, whereas reflection intensities are mostly disregarded, except for the precise localization of reflections and the discrimination of proper reflections from artifacts by considering peak profiles (post-refinement is outside of scope in the present context). If the unit cell parameters are known and accurate positions of the Bragg peaks have been extracted from a diffraction image, the orientation of the crystal lattice can in many cases be calculated within a fraction of a second by using one of the readily available indexing tools.

Since UED experiments typically lack precise information about the position of individual Bragg peaks, Garfield's first design principle (P1) is to use as much of the available information as possible. Hence, the discriminating information in reflection intensities should not be ignored. However, the application of this principle has consequences that go well beyond the inclusion of intensity data, as discussed below. The second principle (P2) requires that the fuzziness of recorded diffraction data (due to imperfections of the sample, for instance) should be taken into account in the process of finding the best interpretation of the data. First, it can provide additional discriminating information, and second, it allows the (mean) orientation of the crystal lattice to be fitted with standard least squares methods, which would not be possible for diffraction of a perfect crystal of infinite size assuming ideal scattering geometry.

NLS fitting of 20 or more parameters, individually or simultaneously, is challenging. Convergence is not guaranteed, and starting values not too far from the solution are required in order to avoid convergence toward a wrong local minimum. Unsupervised fitting is therefore not practical under such conditions. Hence, Garfield's workflow is not designed to work like a black-box which automatically returns a result when a set of input data is entered, but rather provides a workbench that allows the user to work with multiple parameter sets and interactively explore the parameter space until a satisfactory result is found. This is the second fundamental difference between Garfield and other indexing tools. In essence, Garfield displays and manages a table of parameter sets that can be modified and extended by the user. In the following, this table will be referred to as the “main table of settings,” or the main table for short. Garfield provides various tools to create and modify settings, either by manual editing, by least squares optimization of existing parameter sets (GeoFit), or by systematic screening of crystal lattice orientations (GridScan). The goal is to arrive at a setting that optimally accounts for the totality of the recorded diffraction information.

The quality of different parameter sets can be evaluated by SSR values and crystallographic R-factors (after NLS optimization), and by visual comparison of Garfield's predicted diffraction with the observed diffraction. Although not strictly quantitative, the importance of visual comparison cannot be overstated, since it is often the prediction of fine details in the recorded diffraction, not captured by the input list of reduced data, that makes a particular parameter set most convincing. This can be viewed as a special kind of cross-validation and is consistent with the first design principle mentioned above (P1), which would be violated if only the reduced data were considered.

Garfield consists of a collection of scripts written in the R programming language for statistical computing.³⁹ R is a free software available for common computer platforms. Garfield needs installation of R and additional packages available from the Comprehensive R Archive Network (CRAN). Although most of the code is platform-independent, Garfield currently runs only on a Linux system with the GTK2 and Cairo libraries installed. Detailed step-by-step installation instructions are included in the documentation/user manual provided in form of a TiddlyWiki,⁴⁰ the file GarfieldWiki.html, which is part of the Garfield distribution.

A Garfield project is defined by (1) a table of reflection positions and intensities extracted from either a single diffraction image to be indexed, or a tilt series of diffraction images; (2) the crystal structure (in form of a CIF file) which must be known a priori; and (3), optionally, a diffraction image (one or several in the case of a tilt series) in TIFF, JPEG, or PNG format for visual comparison and verification. Here and in the following, the simplest case of a project for indexing a single diffraction image will be assumed. The generalization to multi-image projects, e.g., for tilt-series obtained by rotating the sample about an axis perpendicular the incident beam, is self-explanatory.

Although all calculations are based exclusively on the reduced data (1) and the crystal structure (2), and thus, diffraction images (3) are not strictly necessary for Garfield to work, it is strongly recommended to attach representative diffraction image files to a project (see Sec. II). Full functionality is only attained if a diffraction image has been assigned. For instance, if a diffraction image is assigned, it is also possible to define a mask in order to exclude data from unreliable areas of the detector, such as the shadow of a beamstop. Furthermore, the quality of the results can be easily cross-checked within Garfield by comparing measured with simulated diffraction images. Otherwise, visual comparison is only possible by means of diffraction patterns (in Garfield's parlance), i.e., by two-dimensional scatter plots showing only the positions and integrated intensities of (observed or calculated) Bragg spots.

In Garfield's terminology, which is also adopted in the accompanying documentation and the remainder of this report, a diffraction pattern is a two-dimensional graphical representation of the reduced data (“obs”) or the predicted version of them (“clc”), obtained by plotting circular disks of different size and gray level (to indicate the intensities $Iifoo$⁠) centered at coordinates $(xifoo,yifoo)$⁠, “foo” = “obs” or “clc.” In contrast, a diffraction image corresponds to a matrix of gray values representing the pixel-by-pixel distribution of diffraction intensities, as if recorded with a two-dimensional array detector. Garfield provides options to view the observed diffraction images or patterns and to compare them with simulated counterparts.

The main window of Garfield displays a list of parameter sets (“settings”) for simulating the diffraction data of the current project (see Fig. S1 in the supplementary material).⁷³ If the current project has just been created, the list consists of a single line, representing a dummy parameter set. In this case, in order to start working with the new project, the user must first create a new parameter set by editing the default parameters and assigning them realistic or at least meaningful values. Most of the parameters (including the orientation of the crystal lattice) should be reasonably well known from the experimental setup, but sometimes this is not the case. Other parameters, such as the mosaicity and the average size and shape of the crystallites, are usually not known in advance and must be fitted to the data by using the GeoFit tool.

The main window is used to list and manage a growing collection of parameter sets used to explore the parameter space in order to find the optimal description of the diffraction experiment. New parameter sets can be added by manually editing existing sets or by saving the results of GridScan and GeoFit runs.

A complete set consists of model parameters (Table I) and control parameters (Table II). The former describe the diffraction geometry and sample properties. Many of them can be fitted with GeoFit. Control parameters include binary flags and numerical parameters that define details of how diffraction images are predicted and fitted, such as the type of modeling to be used (“standard” or “ANISO,” see Sec. IV D), the resolution range, or the ambit radius (radius of circular areas around observed diffraction spots used in the definition of the cost function, see Sec. IV A), and whether all predictions or only matching ones (“hits”) should be considered (for further explanations see Sec. IV).

Any setting can be selected as the “active setting” by double-clicking on a row in the main table. Only one setting can be “active” at a time. Most of the options available via drop-down menus can only be used with the active setting defined. The parameters of the active setting are used when simulating diffraction patterns and diffraction images (menu “Display”), when searching for possible crystal orientations (menu “GridScan”), or as starting values for NLS parameter fitting (menu “GeoFit”).

Once a new setting with more realistic parameter values has been created by manual editing, it is advisable to compare the diffraction pattern predicted by this setting with the observed pattern. If the predicted pattern shows resemblance to the observed pattern, it might be possible to proceed directly to the final step of optimizing the parameters by NLS fitting with GeoFit. Often this is not the case, and for the parameter values that are not known precisely enough, better start values for GeoFit must be found. The most critical parameters are clearly the crystal orientation parameters, because even a small misalignment can completely change the diffraction pattern. Other parameters that are typically affected by large uncertainties, such as mosaicity and effective domain size, are less critical, because the diffraction changes continuously with these parameters.⁴¹ 

The GridScan tool is an indispensable aid for finding better candidates for the orientation of the crystal lattice by performing a grid search over the entire $SO(3)$ rotation space.⁴² When the tool is activated, a new window opens for setting up, starting, stopping, and monitoring the progress of grid searches (Fig. 2). In each run, the parameters of the active setting (the one that was active when GridScan was launched), except the orientation parameters and the resolution range, are used to compute a figure of merit (FOM)⁴³ for up to 17 694 720 (⁠ $=49 152×360$⁠) orientations, by comparing the predicted patterns with the “observed” pattern.

In Garfield, the orientation of the crystal is defined by two polar angles, $Θ$ and $Φ$⁠, which define the direction of the incident beam relative to the crystal lattice, and a third angle, $Ψ$⁠, for the rotation of the crystal around the incident beam.⁴⁴ The grid is constructed from 49 152 predefined directions of the beam $(Θ,Φ)$⁠, each one combined with 360 rotations sampling $Ψ$ in steps of 1°. The beam directions and the order in which they are visited are defined according to the principles of HEALPix (Hierarchical Equal Area isoLatitude Pixelation⁴⁵). This sampling results in a most uniform and efficient⁴⁶ coverage of the sphere with an angular resolution of 0.92°. The angular distance of any arbitrary orientation from the nearest grid point is less than 1°.

For each pair $(Θ,Φ)$⁠, represented by a point on the unit sphere, the maximum value of the FOM obtained by scanning over $Ψ$ is plotted in two-dimensional gray-scale maps showing the projections of the upper and lower hemispheres onto the x and y plane.⁴⁷ Good candidate orientations will show up as dark spots in these maps. On a modern multi-core workstation, a full grid scan takes only a few minutes to half an hour, depending on the number of diffraction spots used (typically 50–100). This number depends on the resolution range used for the grid scan, which can be set in the GridScan window before starting a scan. Also, the number N of beam directions to evaluate has to be set before starting a run (⁠ $N≤49152$⁠). A running grid scan can be aborted at any time (and resumed or extended at a later time, if this seems promising). Once a grid scan is finished (or has been terminated), the n best orientations on the FOM ranking can be exported to Garfield's main table, creating n new settings.

After creating a new setting (either by manual editing or using GridScan), it is usually necessary to optimize the current parameters by NLS fitting using the GeoFit tool (for more details, see the supplementary material, Fig. S2).⁷³ The last step in a project should always be to run GeoFit to optimize the parameters and get a summary of the final results.

The fitted values of all model parameters and the control parameters used for fitting can be saved for further use. Saving the results creates a new setting and adds a new entry to the Garfield's main table. In addition, several files are written to the file system for documentation, including (1) a listing of all parameters with their values before and after fitting, (2) a corresponding list of various kinds of SSR values and R factors, and (3) a comprehensive table of data related to Bragg reflections, combining the input coordinates and intensities of observed Bragg spots, the corresponding calculated intensities and centroid positions, and—for each spot—the proposed reflection indices hkl, expected intensities, and excitation errors (distances from the Ewald sphere) of up to five Bragg reflections contributing to the Bragg spot. A pdf file with a graphical representation of the most relevant information contained in the table (3) is also generated (see Fig. 3).

Note that in the example of Fig. 3, none of the observed spots is caused by a single Bragg reflection, but all spots are superpositions of multiple reflections. By considering spots with a dominant and one or several minor contributions, two extreme cases can be distinguished: (i) The structure factors are comparable in magnitude, but the excitation errors of the minor contributions are large compared to the dominant reflection. (ii) The excitation errors are comparable, but the structure factors are very different in magnitude. In case (i), but not in case (ii), it may be justified to neglect the contribution of spurious reflections to the temporal intensity variations measured in UED experiments.

The purpose of the two core tools of Garfield, GridScan and GeoFit, is to find model parameters that reproduce the input data, the positions, and intensities of recorded Bragg spots (the reduced data). Bragg spots are either individual Bragg reflections or superpositions of them. In Garfield's terminology, the term diffraction pattern refers to the graphical representation of the reduced data. Thus, the two terms are essentially equivalent.

Predicting a diffraction pattern for comparison with the observed pattern requires (1) a crystallographic model describing the diffraction by the sample, and suitable algorithms to (2) estimate the expected diffraction and to (3) convert the result into a list of triplets $(xioclc,yioclc,Iioclc)$ analogous to the input list of reduced data $(xioobs,yioobs,Iioobs)$⁠, where index $io=1,… Nobs$ enumerates the observed Bragg spots. In order to be of practical value for NLS fitting of the model parameters, the model must be simple to allow fast computation. Also, the focus is on diffraction patterns (rather than diffraction images) because the algorithms computing the cost function in NLS fitting should be as fast as possible.

The potential loss of precision of the predicted spot positions $(xioclc,yioclc)$ due to various simplifications and approximations in the model is compensated by additional restraints imposed on the model by considering the spot intensities as input data. The goal is to find the solution corresponding to the global minimum of the cost function. Taking intensities into account increases the contrast between different solutions, even if the calculated intensities $Iioclc$ suffer from large uncertainties, but the correlation between estimated and observed intensities is positive. If a well-separated solution can be identified, the orientation of the crystal lattice is usually quite robust to errors in the calculated intensities, so that the assignment of Laue indices is not affected. Thus, the effect of considering intensities can be understood as the application of a filter that eliminates unlikely solutions. The effect of intensity uncertainties is to reduce the selectivity of the filter, rather than to bias the solution toward a completely different set of parameters, i.e., a solution that would produce a wrong indexing.

The following subsections describe how the above requirements (1–3) are implemented in Garfield. We start with (3), explaining how the predicted diffraction pattern is derived from the simulated diffraction. This provides an opportunity to introduce the concept of spot ambits, which is essential for understanding the whole approach and helps to put the limitations of the diffraction model into the right perspective. This first subsection (A) is followed by the definition of the cost function used in NLS fitting (B). The last two subsections describe the model itself (C) and how it is used to simulate the diffraction (D).

In the kinematical theory of diffraction, a perfect crystal in a parallel, monochromatic beam will diffract only in well-defined directions and only if the crystal is oriented so that at least one point of the reciprocal lattice (in addition to the origin) fulfills the Bragg condition by intersecting with the Ewald sphere. If the crystal consists of small domains with slightly different orientations and the beam is not perfectly coherent, RLPs near the Ewald sphere can also give rise to diffraction. This can be viewed as an expansion of the RLPs from points of zero extent to distributions over a finite volume. When treated as distributions in reciprocal space, virtually all RLPs (up to a certain resolution) can contribute to diffraction, but most of them, however, with an intensity that is zero or near-zero.

For now, let us assume that the contributions of all RLPs, enumerated by their Laue indices hkl, have been calculated (as described in Secs. IV C and IV D) and are available as a table of centroid coordinates, integrated intensities, and corresponding Laue indices: $xjp, yjp, Ijp$⁠, and $(hkl)jp$ with $jp=1,… Np$⁠. The number $Np$ of predicted reflections is usually much larger than the number of observed spots $Nobs$⁠.

For quick comparison with the reduced data (the input diffraction pattern), Garfield provides the option to calculate the values in Eqs. (5) and (6) for any setting in the main table and to display the result as simulated diffraction pattern.

The parameters of a selected setting can be fitted with the GeoFit tool by NLS minimization of a cost function S, which is defined as the sum of squared differences between calculated and observed quantities that are related to the intensities and coordinates of diffraction spots [cf. Eq. (1)]. The exact form of the cost function can be adjusted by using a set of control parameters. One of the control parameters is $wP$ already introduced in Eq. (1). Other control parameters are the minimum and maximum d-spacing of the resolution range to be considered.⁵⁰ 

The $wio$ in Eq. (8) are individual weights, which in principle should be read from the input list of reduced data. However, weights (or sigma values) assigned to observed intensities are not taken into account in the present version of Garfield. Instead, the weights are determined as $wio=(Iioobs)−12$⁠, which has the effect to weighing down high intensities relative to low intensities, similar to the effect of setting $target=“sqrt(I)”$⁠.

It should be noted that position information present in the reduced data enters the cost function via the ambits, even if $wP=0$ (the default) and, thus, $S(P)=0$⁠.

All simulations in Garfield are calculated within the kinematical theory of diffraction. Multiple scattering, which can pose a serious problem for structure determination by electron diffraction, is ignored. For nearly perfect crystals, the dynamical theory of diffraction can in principle be used to cope with the effect of multiple scattering. This is a much more difficult (though not impossible⁵²) task for typical samples used in transmission UED, where the sample quality is less perfect due to the reasons described in the Introduction. Thus, the application of dynamical theory in Garfield is not practical due to the enormous computational effort that would be required, and it is questionable whether the results would improve significantly.

Fortunately, the method used in Garfield for indexing UED patterns can tolerate quite large errors in the estimated diffraction intensities, as explained at the beginning of this section (Sec. IV). Furthermore, there are solid reasons to assume that dynamical effects are less pronounced in highly textured crystal, i.e., when the coherently diffracting domains are small and the mosaic spread is large^53,54 as is often the case due to the typical limitations of UED experiments. This has recently been confirmed by simulations and comparison of electron diffraction calculated in kinematical and dynamical theory.⁵² According to the numerical results for single-crystal gold, application of the dynamical diffraction theory was indispensable for the accurate determination of the lattice temperature. However, even in this extreme case (closely packed heavy atoms) the errors of the kinematical intensities, measured by the R factor between kinematical and dynamical simulation, did not exceed much the 30% level, except for rather small tilt spreads (⁠ $σθ<20 mrad$⁠). An R value of 30% indicates that kinematical intensities are subject to large errors, but they nevertheless contain useful information that Garfield will exploit to determine the correct indexing.

Ignoring dynamical diffraction effects reduces the description of electron diffraction to a mere problem of considering the relative phases of waves scattered by the distribution of matter in the crystal, in strict analogy to x-ray diffraction. This leads to the Bragg or Laue conditions of diffraction, which can be visualized with the Ewald sphere construction. The transition from x-ray to electron diffraction is then essentially a matter of replacing x-ray atomic scattering factors with electron scattering factors and taking into account that for electrons, the radius of the Ewald sphere, $1/λ$⁠, is much larger than for X-rays and almost “plane-like.” Garfield uses the parameterizations of electron scattering factors of atoms and ions by Peng et al.^55,56

The restriction to very thin samples relaxes the diffraction conditions to some extent, so that points of the reciprocal lattice can give rise to reflections even if the Laue conditions are not exactly fulfilled, provided that the excitation error, the distance of an RLP from the Ewald sphere, is not too large. Usually, this situation is interpreted in a different way: Instead of relaxing the Laue conditions, the RLPs are replaced by needle-shaped distributions parallel to the surface normal of the crystal sheet, and diffraction occurs, if a “reciprocal lattice rod” (relrod) intersects with the Ewald sphere. This is consistent with the idea, that the reciprocal lattice weighted with the structure factor represents the Fourier transform of a perfect crystal of infinite size, $F∞$⁠, whereas the Fourier transform of a finite crystal is the convolution of $F∞$ with the Fourier transform of the crystal's shape function. For a plane-parallel plate of thickness t perpendicular to the incident beam, this leads to relrod profiles oscillating like a $sinc$ function. The intensity measured for a particular RLP with indices hkl varies as $sinc2(π t Δh)$ with the distance $Δh$ of the RLP from the Ewald sphere (measured in units of $h=1/dhkl$⁠, the inverse of the d-spacing).

In UED experiments, the finite size of the crystal is not the only, and often not the most important reason why non-zero diffraction intensities occur for RLPs even if they are nominally not located on the Ewald sphere. Because of the large area of the diffracting volume (compared to its thickness) and the non-perfect quality of most samples, the orientation of the crystal lattice is not exactly the same in different parts of the sample. Thus, UED samples should rather be viewed as an ensemble of crystalline domains with slightly different orientations distributed around the average orientation of the crystal lattice. This is similar to the orientation distribution of mosaic blocks introduced by Darwin⁵⁷ to quantitatively explain x-ray diffraction of single crystals. Although the orientation distribution in UED experiments can be much wider (typically several degrees) than the typical mosaic spread of good single crystals (usually $≪0.1°$⁠), the effect of orientation disorder will be referred to as mosaicity in the following. Other effects that contribute to the appearance of reflection spots in UED images are the divergence of the incident beam and, to a lesser extent, its finite energy bandwidth. In Garfield, all these effects are accounted for by following the example of the relrods introduced to explain the finite crystal thickness effect: RLPs (represented by Dirac delta functions) are replaced by smooth distributions, which should lead to correct results, at least approximately, if the Ewald sphere construction is used without modification.

To allow simple and fast computation of the combined effect, Garfield assumes normal distributions to describe the four effects. Gaussian approximations have been used in various formulations in serial snapshot crystallography to estimate reflection intensities and partialities (see Brehm et al.,⁵⁸ and the references cited therein). In Garfield, the simplest case is assumed: each individual effect is described by a single Gaussian distribution, and the distributions are mutually independent. Then, the combination of any two effects corresponds to the convolution of two Gaussians, which in turn is a Gaussian distribution. Furthermore, the marginal and conditional distributions of multivariate Gaussians are also normal distributions. This allows a simple analytical formula to be used in many cases, while more realistic distributions would require time-consuming numerical integration. A few remarks are in order for all four effects.

Given the fact that the intensity profile of the relrods for a thin plate [Fig. 4(b)] oscillates like the square of the $sinc$ function, using normal distributions to describe the finite size effect may seem unphysical. However, the set of coherently diffracting domains (crystallites) of a typical UED sample varies not only in orientation but also in lateral extent and thickness. Averaging over crystallites of different thicknesses will fill in the minima of the intensity profile, $PI(z)$⁠, and results in smoothing, making the profile sufficiently normal. The smoothing will be even more pronounced if limitation of the domains in lateral directions (perpendicular to the beam) cannot be neglected. In this case, the relrods must be described by three-dimensional distributions of finite extension perpendicular to their main direction. An extreme case of diffraction by isometric particles is represented in Fig. 4(d) by the intensity profile for the solid sphere.

The intensity profile of the RLPs due to the finite size effect is described by a general three-dimensional normal distribution with six independent parameters: three standard deviations $σishp$ along the principal axes (⁠ $i=1,2,3$⁠) and three Euler angles. Degenerate distributions (one or two of the $σishp$ equal to zero) are allowed. Garfield is able to handle the case of a thin plate, all other broadening effects being negligible. However, due to the restriction to Gaussian profiles, the NLS fitting of thickness t is biased toward smaller values than the true average $t¯$⁠, since in this way the slow decay of the side maxima of the $sinc$ function can be compensated by increasing the width of the Gaussian profile, $σq∝t−1$⁠. In other words, Garfield cannot be used to fit the geometric thickness of the sample because the parameter t is correlated with, but not identical to, an effective thickness.

In contrast to the finite size effect, which is the same for all RLPs, the broadening of RLPs due to the orientation distribution of mosaic blocks⁵⁹ increases linearly with the resolution, $q=1/d$⁠. Here, d is the d-spacing, $q=2 sinθ /λ$⁠, $θ$ the Bragg angle, and $λ$ the de Broglie wavelength of the electrons. Often the mosaicity can be quantified by a single angle, $σω$⁠, the standard deviation of rotation angles $ω$ between the mean orientation and the actual orientations of the mosaic blocks (or “domains”). This assumes that the mosaic spread is isotropic.

In typical UED experiments, a significant contribution to the mosaic spread is due to flatness imperfections resulting from the mechanical instability of thin plates and difficulties in sample preparation and mounting. Orientation distributions owing to such imperfections are not necessarily isotropic. For example, uniform bending of a flat plate would result in a systematic variation in local orientation (represented by the tangent planes), which is manifestly anisotropic.

Garfield accounts for the anisotropy of the mosaic spread by assuming that the rotation vectors $ω→$⁠, which represent the rotations of the domains from the mean orientation (averaged over the whole sample) to the actual orientation, are distributed according to a three-dimensional normal distribution.^60,61 In general, this orientation distribution, $Pω→$⁠, is described by six independent parameters, e.g., the standard deviations $σimos (i=1,2,3)$ along the principal axes and three Euler angles defining the orientation of the principal axes relative to the laboratory coordinate system (orthogonal axes X, Y, and Z with Z parallel to the incident beam).

The distribution $Ph→$ of the domains' reciprocal lattice vectors $h→$ induced by the rotations⁶² is determined by $Pω→$⁠. Viewed as a three-dimensional distribution in reciprocal space, $Ph→$ is degenerate, as it occupies a two-dimensional surface given by a sphere of radius $h=|h→|$ centered at the origin. Other effects on the intensity profile set aside, diffraction conditions require to consider the intersection of this sphere with the Ewald sphere. Thus, diffraction intensity contributions due to mosaicity will appear as circular arcs. The intensity distribution along the arcs is defined by the conditional distribution obtained by restricting $Ph→$ to the (curved) line of intersection. From this, the expected intensity profiles on the detector can be derived, and their total intensities (integrated along the arcs) as well as the centroid positions can be calculated.

The exact calculation would be cumbersome and would require numerical integration of functions in three dimensions for all RLPs $h→=i a→⋆+j b→⋆+k c→⋆$ of interest (⁠ $a→⋆$⁠, $b→⋆$⁠, and $c→⋆$ basis vectors of the reciprocal lattice). Garfield exploits the fact that $Pω→$ is a normal distribution (by assumption) and introduces an approximation which is only valid if all relevant rotation angles are sufficiently small. In the following, the approximation is first motivated by considering the case where all rotations are infinitesimally small; subsequently, the question is addressed as to how large the angles of rotation may be so that the approximation can still be used for practical purposes.

The first step in calculating the expected intensity profile of the reflection corresponding to a given RLP $h→$ is to transform the density function of $Pω→$ from the laboratory coordinate system to a local Cartesian coordinate system $(ω1,ω2,ω3)$ with the third axis parallel to $h→$⁠. The other two axes are chosen to be perpendicular (⁠ $ω1$⁠) and parallel (⁠ $ω2$⁠) to the plane, which is spanned by $h→$ and the wavevector $ki→∥Z$ of the incident beam (⁠ $|ki→|=1/λ$⁠), see Fig. 5. The covariance matrix of the transformed density function will be denoted by $Σω→(123)$⁠.

From the marginal distribution $Pω→⊥⊥h→$⁠, a density function for the rotated $h→$ vectors can be derived straightforwardly. Within the infinitesimal rotations approximation, the distributions $Ph→$ and $Pω→⊥⊥h→$ are interchangeable, and the density function in Eq. (12) is also valid for $h→$ if variables $(ω1,ω2)$ are interpreted as a special parameterization of $h→$⁠.

In the tangent plane approximation, $ρ1$ and $ρ2$ can be identified with $h1$ and $h2$⁠, the first two coordinates of $h→$ in the rotated coordinate system.

The expected intensity profile of the reflection corresponding to the RLP $h→$ can be obtained (up to a scaling factor containing the squared amplitude of the structure factor) by calculating $Ph→(h1,h2)$ along the trace of the Ewald sphere in the tangent plane. To a very good approximation, this trace is a straight line. The result is a normal distribution whose mean value (center position) and variance can be calculated from $Σh→$ by using another transformation of the coordinate system, such that the trace of the Ewald sphere is parallel to one of the new axes.

For non-infinitesimal rotations, the tangent plane approximation breaks down, and the curvature of the support of $Ph→$⁠, i.e., the sphere of radius h around the origin, has to be taken into account. This is done in Garfield by applying geometric corrections that should perform well for all practical purpose.⁶⁴ Hence, the most critical step in estimating diffraction intensities is the infinitesimal rotations approximation used for the transition from the original distribution $Pω→$ to its marginal distribution $Pω→⊥⊥h→$ and the identification of $Pω→⊥⊥h→$ with $Ph→$⁠.

The limitations of this approximation are illustrated in Fig. 6. The density of the marginal distribution at any point of the $(ω1,ω2)$ plane is obtained by integration of $Pω→$ along the $ω3$ axis, which is parallel to a certain RPL's $h→$ vector. However, this is not exactly what is required as rotation vectors $(ω1,ω2,0)$ and $(ω1,ω2,ω3≠0)$ do not move the RLP exactly to the same position, unless $ω3$ is infinitesimal. In order to integrate densities of rotations that move the RLP to identical positions, $Pω→$ should be integrated along curves such that all points of these curves correspond to rotations that result in the same position of the RLP. Projections of such curves to the $(ω1,ω2)$ plane are shown in Fig. 6.

One observation that follows directly from looking at Fig. 6 is that the marginal distribution intermixes orientations (⁠ $ω→$ vectors) within angular ranges that increase almost linearly with $ω3$⁠, leading to a certain alteration of the correct distribution. The figure also shows that intermixing of nonequivalent orientations has virtually no effect if the original distribution is isotropic, and that errors due to intermixing increase with the degree of anisotropy. Except in extreme cases of anisotropy, reducing the marginal distribution to a line in the $(ω1,ω2)$ plane (for example the line through the black dots in Fig. 6), standard deviations $σimos$ up to 5° (i.e., 3 sigma < 15°) should be acceptable, considering that Garfield's requirements for the precision of position and intensity estimates are not very high. Even sigma values in excess of 10° may be tolerable if the anisotropy is moderate.

In Garfield, the angular spread of the incident wavevectors, $k→i$⁠, is modeled as a rotationally symmetric normal distribution with sigma values $σx=σy=σdivg$⁠. The effect of a finite beam divergence on the diffraction image is treated approximatively by splitting the effect in two parts and treating them separately. (i) The first effect is that more reflections can appear because more RLPs may “touch” the bundle of Ewald spheres corresponding to the distribution of $k→i$ vectors, and the intensities of all reflections change according to the change of “mean” or “effective” excitation errors. (ii) The second part concerns the reflections' positions and profiles in the detector plane, which are affected by beam divergence due to the range of projection directions associated with the angular spread of the $k→i$ vectors.

The first effect is taken into account by replacing the distribution $Pω→$ of orientation vectors (the “mosaicity”) by the convolution of $Pω→$ with the distribution describing the divergence of the $k→i$ vectors.⁶⁵ Thus, this part is treated by replacing the mosaicity with an “effective mosaicity,” while the Ewald sphere construction remains unchanged.

In the calculation of the cost function (for parameter fitting), the second effect is ignored. The cost function only depends on the center positions and integrated intensities of predicted reflections. The profiles are irrelevant, and the shift in center positions is negligible if the beam divergence is significantly smaller than the mosaicity, which is the case in typical UED experiments. Nevertheless, for simulating the diffraction image, the broadening of reflection profiles by $σdivg$ is approximately taken into account.

An energy distribution around the mean electron energy leads to a distribution of $ki=1/λ$ values and thus to a distribution of Ewald spheres with radii $ki$⁠. Again, Garfield uses a normal distribution (sigma value $σbwdth$⁠) to describe the effect of a finite bandwidth. Similar to beam divergence, the energy bandwidth affects (i) which RLPs give rise to a reflection and with what intensity, and (ii) where in the diffraction image the reflections appear and what their intensity profile is. Garfield treats the effect of energy distribution in an approximative manner, by separating the two aspects. This leads to another increase in the effective mosaicity by convolution with the distribution describing the energy spread. In this case, however, the modification of $Pω→$ depends on the resolution, i.e., the distance of the RLP from the origin, $|h→|$⁠. For the same reason as in the treatment of beam divergence, the second effect (ii) is considered only when simulating diffraction images (where reflection profiles should be as realistic as possible). For parameter fitting and predicting diffraction patterns, the second effect is ignored.⁶⁶ 

Using normal distributions to describe the four effects discussed in Sec. IV C has the advantage that predicting reflection intensities and positions is largely reduced to algebraic manipulations of covariance matrices. The combination of beam divergence and energy spread with the orientation distribution of coherently diffracting domains is not a major complication and can be achieved by convolution, resulting in an effective mosaicity. This is unproblematic because it describes the broadening of the RLPs' diffracting profiles due to three effects that are all of the same kind: All three are related to relaxing one of the conditions of an ideal diffraction experiment with scattering geometry defined by a unique direction of the incident beam, a unique energy of the electrons, and a unique orientation of the crystal lattice.

The finite size effect is of a different kind. Unless the (average or effective) shape function is spherically symmetric, it is not sufficient—in principle—to know $Ph→(h→)$⁠, the density of RLPs at any point $h→$⁠, without also considering the orientation of the corresponding subset of crystallites, which deviates from the average orientation of the crystal lattice by virtue of $ω1=−ρ2/h$⁠, $ω2=ρ1/h$⁠, and $ω3$⁠. Within the tangent plane approximation, rotation of the shape function (represented by normal distribution $Pshp$⁠) due to $ω1$ and $ω2$ can be neglected. Thus, for this part, a simple convolution of $Ph→$ with $Pshp$ would suffice. Extending this approach to larger values of $ω1$ and $ω2$ should be possible by means of the geometric corrections mentioned in Sec. IV C 2.⁶⁷ 

With regard to $ω3$⁠, an exact treatment is computationally more expensive, since all information about $ω3$ rotations is lost when calculating $Ph→(ω1,ω2)$⁠, the marginal distribution of $Pω→(ω1,ω2,ω3)$⁠. The exact treatment would involve three steps: instead of using $Ph→(ω1,ω2)$⁠, consider the conditional distributions $Pω→(ω1,ω2 | ω3) dω3$⁠, calculate their convolution with $Pshp$ rotated by $ω3$⁠, and numerically integrate the results over $ω3$⁠. We considered such a protocol, but the results did not justify the extra computational costs associated with it, and in the current version of Garfield, the rotation of the shape function due to $ω3$ is ignored.⁶⁸ 

Convolution of the bivariate (degenerate trivariate) distribution $Ph→$ with the shape function results in a non-degenerate trivariate distribution. This poses no further problems. The result is a two-dimensional intensity profile with non-zero width in radial direction. The radial width is only determined by the shape function.

Anisotropy of the orientation distribution appears to be a common phenomenon in UED, but it is not always a dominant feature. Whenever possible, the general treatment described above should be avoided by using a simpler approach that permits much faster calculations by ignoring mosaic anisotropy. In the GeoFit tool and for simulating and displaying diffraction patterns, Garfield offers the option to choose between two methods for predicting position and intensity of reflections: the method presented in Sec. IV C and the previous paragraphs (“ANISO modeling”), and a simplified version called “standard modeling.”⁶⁹ 

In standard modeling, the four factors that determine the RLPs' diffracting power profiles are all described by Gaussian functions with standard deviations $σshp$⁠, $σmos$⁠, $σdivg$⁠, and $σbwdth$⁠. Here, the orientation spread measured by $σmos$ is assumed to be the same in all directions (i.e., isotropic), while the shape profile measured by $σshp$ is restricted to a certain direction (“relrods”). Standard modeling assumes that the coherently diffracting domains of the sample can be viewed as thin plates. Their finite extent in lateral directions is ignored.⁷⁰ 

Finally, $σ*$ is used to estimate the integrated intensity and centroid position of the reflection associated with the RLP. Since $σ*$ involves two fundamentally different effects, this cannot be done in a straightforward manner. Consider the two extreme cases: (1) $σshp≪h σMOS$ and (2) $σshp≫h σMOS$⁠, as illustrated in Fig. 7.

1. If broadening of the RLP profile due to the shape function can be neglected, the profile is determined by the (effective) orientation distribution, so $h→$ is spread over the surface of the sphere $Sh$ with radius h around the origin. Without introducing further assumptions, the best guess for the centroid position of the expected reflection corresponds to the point $A$ on the intersection of the Ewald sphere with $Sh$ and the plane that is spanned by $k→i$ and $h→$⁠. Let $eMOS$ be the distance from the RLP to $A$⁠, measured along $Sh$ in this plane. Then the intensity can be estimated by the value of the profile at $eMOS$ (up to a factor containing the structure factor amplitude squared).

2. If the effective mosaicity is negligible compared to the finite size effect and the shape function is that of a thin plate, $σ*$ represents the diffracting power profile of relrods, which are usually more or less parallel to the (mean) surface normal of the crystal foil. The expected centroid position of the reflection is then determined by the point $B$ where the straight line through the tip of $h→$ parallel to the direction of the relrods hits the Ewald sphere. The estimated intensity is proportional to the value of the Gaussian density profile at $eshp$⁠, the distance between $B$ and the tip of $h→$⁠.