Automated Crystal Orientation Mapping in py4DSTEM using Sparse Correlation Matching

Colin OphusSteven E ZeltmannAlexandra BruefachAlexander RakowskiBenjamin H SavitzkyAndrew M MinorMC Scott

Crystalline materials used in technological applications are often complex assemblies composed of multiple phases and differently oriented grains. Robust identification of the phases and orientation relationships from these samples is crucial, but the information extracted from the diffraction condition probed by an electron beam is often incomplete. We therefore have developed an automated crystal orientation mapping (ACOM) procedure which uses a converged electron probe to collect diffraction patterns from multiple locations across a complex sample. We provide an algorithm to determine the orientation of each diffraction pattern based on a fast sparse correlation method. We test the speed and accuracy of our method by indexing diffraction patterns generated using both kinematical and dynamical simulations. We have also measured orientation maps from an experimental dataset consisting of a complex polycrystalline twisted helical AuAgPd nanowire. From these maps we identify twin planes between adjacent grains, which may be responsible for the twisted helical structure. All of our methods are made freely available as open source code, including tutorials which can be adapted to perform ACOM measurements on diffraction pattern datasets.

Keywords:automated crystal orientation mapping (ACOM)four-dimensional scanning transmission electron microscopy (4D-STEM)nanobeam electron diffraction (NBED)open-source softwarescanning electron nanodiffraction (SEND)


Polycrystalline materials are ubiquitous in technological applications. An ideal crystal structure can be fully defined with a small number of parameters: the 3 vectors defining its unit cell, and the position and species of each atom inside the unit cell . To fully describe crystalline materials in the real world however, we require a description of both the crystal lattice, and all defects present in a given material. These include point defects such as dopants, vacancies, or interstitials , line defects such as dislocations , planar defects including internal boundaries and surfaces , and finally volume defects such as local strain fields . One large subset of crystalline materials are polycrystalline phases, which consist of many small crystalline grains, arranged in either a random or organized fashion. Many material properties such as mechanical strength , optical response , or thermal or electrical conductivity are strongly modulated by the density and orientation of the boundaries between crystalline grains . Thus characterizing the orientation of polycrystalline grains is essential to understanding these materials.

The two primary tools used to study the orientation of polycrystalline materials are electron backscatter diffraction (EBSD) in scanning electron microscopy (SEM), and transmission electron microscopy (TEM). EBSD can measure the orientation of crystalline grains with very high accuracy, but has limited resolution and is primarily sensitive to the surface of materials . Alternatively, we can directly measure the atomic-scale structure and therefore the orientation of polycrystalline grains, either by using plane wave imaging in TEM , or by focusing the probe down to sub-atomic dimensions and scanning over the sample surface in scanning TEM (STEM) . This is possible due to the widespread deployment of aberration correction for both TEM and STEM instruments . Atomic resolution imaging, however, strictly limits the achievable field-of-view, and requires relatively thin samples, and thus is primarily suited for measuring polycrystalline grain orientations of 2D materials .

Another approach to orientation mapping in TEM is to use diffraction space measurements. For crystalline materials, diffraction patterns will contain Bragg spots with spacing inversely proportional to the spacing of atomic planes which are approximately perpendicular to the beam direction (described by both the Laue condition and Bragg equations ). To generate a spatially-resolved orientation map, we can focus a STEM probe down to dimensions of 0.5 to 50 nm, scan it over the sample surface, and record the diffraction pattern for each probe position, a technique referred to as nanobeam electron diffraction (NBED) , scanning electron nanobeam diffraction (SEND) , or four dimensional-scanning transmission electron microscopy (4D-STEM) (we choose this nomenclature for this text) due to the 4D shape of the collected data . 4D-STEM experiments are increasingly enabled by fast direct electron detectors, as these cameras allow for much faster recording and much larger fields of view .

By performing template matching of diffraction pattern libraries on 4D-STEM datasets, we can map the orientation of all crystalline grains with sufficient diffraction signal. This method is usually named automated crystal orientation mapping (ACOM), and has been used by many authors in materials science studies . ACOM experiments in 4D-STEM are highly flexible; two recent examples include implementing ACOM measurements in liquid cell experiments, and adapting the ACOM method to a scanning confocal electron diffraction (SCED) experimental configuration. ACOM is also routinely combined with precession electron diffraction, where the STEM beam is continually rotated around a cone incident onto the sample, in order to excite more diffraction spots and thus produce more interpretable diffraction patterns . Recently, have combined simulations with machine learning segmentation to map orientations of 2D materials, and have used machine learning methods to improve the resolution and sensitivity of orientation maps by training on simulated data.

In this study, we introduce a new sparse correlation framework for fast calculation of orientation maps from 4D-STEM datasets. Our method is based on template matching of diffraction patterns along only the populated radial bands of a reference crystal's reciprocal lattice, and uses direct sampling of the first two Euler angles, and a fast Fourier transform correlation step to solve for the final Euler angle. We test our method on both kinematical calculations, and simulated diffraction experiments incorporating dynamical diffraction. Finally, we generate orientation maps of polycrystalline AuAgPd helically twisted nanowires, and use clustering to segment the polycrystalline structure, and map the shared (111) twin planes of adjacent grains.



#Structure Factor Calculations

The structure factors of a given crystalline material are defined as the complex coefficients of the Fourier transform of an infinite crystal . We require these coefficients in order to simulate kinematical diffraction patterns, and thus we briefly outline their calculation procedure here.

First, we define the reference crystal structure. This structure consists of two components, the first being its unit cell defined by its lattice vectors a\bm{a}, b\bm{b}, and c\bm{c} composed of positions in r=(x,y,z)\bm{r}=(x,y,z), the 3D real space coordinate system. The second component of a crystal structure is an array with dimensions [N,4][N, 4] containing the fractional atomic positions pn=(pa,pb,pc)n\bm{p}_n = (p_{\bm{a}}, p_{\bm{b}}, p_{\bm{c}})_n and atomic number ZnZ_n, for the nnth index of NN total atoms in the unit cell. Together these positions and atomic numbers are referred to as the atomic basis. Because pn\bm{p}_n is given in terms of the lattice vectors, all fractional positions have values inside the range [0,1)[0, 1). The unit cell and real space Cartesian coordinates of the fcc Au structure are plotted in a.

All subsequent calculations are performed in reciprocal space (also known as Fourier space or diffraction space). Thus the next step is to compute the reciprocal lattice vectors, defined by

a=b×ca[b×c]=b×cVb=c×ab[c×a]=c×aVc=a×bc[a×b]=a×bV,\begin{align*} \bm{a}^* &=& \frac{\bm{b} \times \bm{c}}{ \bm{a} \cdot [\bm{b} \times \bm{c}]} = \frac{\bm{b} \times \bm{c}}{V} \nonumber \\ \bm{b}^* &=& \frac{\bm{c} \times \bm{a}}{ \bm{b} \cdot [\bm{c} \times \bm{a}]} = \frac{\bm{c} \times \bm{a}}{V} \nonumber \\ \bm{c}^* &=& \frac{\bm{a} \times \bm{b}}{ \bm{c} \cdot [\bm{a} \times \bm{b}]} = \frac{\bm{a} \times \bm{b}}{V}, \end{align*}

where ×\times represents the vector cross product and VV is the cell volume in real space. Note that this definition does not include factors of 2π2\pi, and therefore all reciprocal coordinates have spatial frequency units.

Next, we calculate the position of all reciprocal lattice points required for our kinematical diffraction calculation, given by

ghkl=ha+kb+lc,\bm{g}_{hkl} = h \, \bm{a}^* + k \, \bm{b}^* + l \, \bm{c}^*,

where hh, kk, and ll are integers representing the reciprocal lattice index points corresponding to the Miller indices (h,k,l)(h,k,l). We include only points where qhkl<qmax|\bm{q}_{hkl}| < q_{\rm{max}}, where q=(qx,qy,qz)\bm{q}=(q_x,q_y,q_z) are the 3D coordinates in reciprocal space, i.e.\ those which fall inside a sphere given by the maximum scattering vector qmaxq_{\rm{max}}. To find all reciprocal lattice coordinates, we first determine the shortest vector given by linear combinations of (a,b,c)(\bm{a}^*,\bm{b}^*,\bm{c}^*), and divide qmaxq_{\rm{max}} by this vector length to give the range for (h,k,l)(h,k,l). We then tile (h,k,l)(h,k,l) in both the positive and negative directions up to this value, and then remove all points with vector lengths larger than qmaxq_{\rm{max}}.

The reciprocal lattice defined above represents all possible coordinates where the structure factor coefficients Vg(q)V_g(\bm{q}) could be non-zero. The structure factor coefficients depend only the atomic basis and are given by

Fhkl=n=1Nfn(ghkl)exp[2πi(h,k,l)pn],F_{hkl} = \sum_{n=1}^N f_n(|\bm{g}_{hkl}|) \exp\left[ -2 \pi \ii (h,k,l) \cdot \bm{p}_n \right],

where fnf_n are the the single-atom scattering factors for the nnth atom, which describe the scattering amplitude for a single atom isolated in space. There are multiple ways to parameterize fnf_n, but here we have chosen to use the factors defined by which are implemented in py4DSTEM\pyFDSTEM. b shows the atomic scattering factor for an Au atom.

We have now defined all structure factor coefficients for a perfect infinite crystal as

Vg(q)={Fhklif q=ghkl0otherwise. V_g(\bm{q}) = \begin{cases} F_{hkl} & \text{if $\bm{q} = \bm{g}_{hkl}$} \\ 0 & \text{otherwise}. \end{cases}

c shows the structure factors of fcc Au, where the marker size denotes the intensity (magnitude squared) of the FhklF_{hkl} values.

#Calculation of Kinematical Diffraction Patterns

Here we briefly review the theory of kinematical diffraction of finite crystals, following . We can fully describe an electron plane wave by its wavevector k\bm{k}, which points in the direction of the electron beam and has a length given by k=1/λ|\bm{k}| = 1/\lambda, where λ\lambda is the (relativistically-corrected) electron wavelength. Bragg diffraction of the electron wave along a direction k\bm{k}' occurs when electrons scatter from equally spaced planes in the crystal, described in reciprocal space as

k=k+ghkl.\bm{k}' = \bm{k} + \bm{g}_{hkl}.

For elastic scattering, k\bm{k}' has the same length as k\bm{k}, and so scattering can only occur along the spherical surface known as the Ewald sphere construction . This expression will almost never be satisfied by a perfect infinite crystal. However, real samples have finite dimensions, and thus the Fourier transform of their lattice will include a shape factor D(q)D(\bm{q}) convolved with each reciprocal lattice point. Thus diffraction can still occur, as long as is approximately satisfied.

If the sample foil is tilted an angle α\alpha away from the beam direction, the vector between a reciprocal lattice point g\bm{g} and its closest point on the Ewald sphere has a length equal to

sg=g(2k+g)2k+gcos(α).s_{\bm{g}} = \frac{-\bm{g} \cdot(2 \bm{k} + \bm{g}) } {2 | \bm{k} + \bm{g} | \cos(\alpha) }.

The sgs_{\bm{g}} term is known as the excitation error of a given reciprocal lattice point g\bm{g}. When the excitation error sg=0s_{\bm{g}}=0, the Bragg condition is exactly satisfied. When the length of sgs_{\bm{g}} is on the same scale as the extent of the shape factor, the Bragg condition is approximately satisfied.

A typical TEM sample can be approximately described as a slab or foil which is infinite in two dimensions, and with some thickness tt along the normal direction. The shape function of such a sample is equal to

D(qz)=sin(πqzt)πqz.D(q_z) = \frac{\sin(\pi q_z t)}{\pi q_z}.

Because this expression is convolved with each reciprocal lattice point, we can replace qzq_z with the distance between the Ewald sphere and the reciprocal lattice point. For the orientation mapping application considered in this paper, we assume that α=0\alpha = 0, and that the sample thickness tt is unknown. Instead, we replace with the approximation

D(qz)=exp(qz22σ2),D(q_z) = \exp\left( -\frac{{q_z}^2}{2 \sigma^2} \right),

where σ\sigma represents the excitation error tolerance for a given diffraction spot to be included. We chose this expression for the shape function because it decreases monotonically with increasing distance between the diffraction spot and the Ewald sphere qzq_z, and produces smooth output correlograms.

To calculate a kinematic diffraction pattern for a given orientation w\bm{w}, we loop through all reciprocal lattice points and use to calculate the excitation errors. The intensity of each diffraction spot is given by the intensity of the structure factor Fhkl2|F_{hkl}|^2, reduced by a factor defined by either or . We define the position of the diffraction spots in the imaging plane by finding two vectors perpendicular to the beam direction, and projecting the diffraction vectors qq into this plane. The result is the intensity of each spot ImI_m, and its two spatial coordinates (qmx,qmy)(q_{m_x},q_{m_y}), or alternatively their polar coordinates qm=qmx2+qmy2q_m = \sqrt{{q_{m_x}}^2 + {q_{m_y}}^2} and γm=\gamma_m = arctan2(qmy,qmx)(q_{m_y}, q_{m_x}). Note that the in-plane rotation angle is arbitrarily defined for kinematical calculations in the forward direction. The resulting diffraction patterns are defined by the list of MM Bragg peaks (qmx,qmy,Im)(q_{m_x},q_{m_y},I_m) or (qm,γm,Im)(q_m,\gamma_m ,I_m).

e shows diffraction patterns for fcc Au, along five different zone axes (orientation directions). Each pattern includes Bragg spots out to a maximum scattering angle of qmax=1.5A˚1q_{\rm{max}} = 1.5 \, \rm{\AA}^{-1}, and each spot is labeled by the (hkl)(hkl) indices. The marker size shown for each spot scales with the amplitude of each spot's structure factor, decreased by using σ=0.02A˚1\sigma = 0.02 \, \rm{\AA}^{-1}.

#Generation of an Orientation Plan

The problem we are solving is to identify the relative orientation between a given diffraction pattern measurement and a parent reference crystal. This orientation can be uniquely defined by a [3×3][3 \times 3]-size matrix m\tensor{\bm{m}}, which rotates vectors d0\bm{d}_0 in the sample coordinate system to vectors d\bm{d} in the parent crystal coordinate system

[dxdydz]=[uxvxwxuyvywyuzvzwz][d0xd0yd0z]d=m d0,\begin{align*} \begin{bmatrix} d_x \\ d_y \\ d_z \end{bmatrix} &=& \left[ \begin{array}{ccc} u_x & v_x & w_x \\ u_y & v_y & w_y \\ u_z & v_z & w_z \\ \end{array} \right] \begin{bmatrix} d_{0_x} \\ d_{0_y} \\ d_{0_z} \end{bmatrix} \nonumber \\ \bm{d} &=& \tensor{\bm{m}} \ \bm{d}_0, \end{align*}

where the first two columns of m\tensor{\bm{m}} given by u\bm{u} and v\bm{v} represent the orientation of the in-plane xx and yy axis directions of the parent crystal coordinate system, respectively, and the third column w\bm{w} defines the zone axis or out-plane-direction. The orientation matrix can be defined in many different ways, but we have chosen to use a ZXZZ-X-Z Euler angle scheme , defined as

m=[C1S10S1C10001][1000C2S20S2C2][C3S30S3C30001],\tensor{\bm{m}} = \left[ \begin{array}{ccc} C_1 & -S_1 & 0 \\ S_1 & C_1 & 0 \\ 0 & 0 & 1 \\ \end{array} \right] \left[ \begin{array}{ccc} 1 & 0 & 0 \\ 0 & C_2 & S_2 \\ 0 & -S_2 & C_2 \\ \end{array} \right] \left[ \begin{array}{ccc} C_3 & -S_3 & 0 \\ S_3 & C_3 & 0 \\ 0 & 0 & 1 \\ \end{array} \right],

where C1=cos(ϕ1)C_1 = \cos(\phi_1), S1=sin(ϕ1)S_1 = \sin(\phi_1), C2=cos(θ2)C_2 = \cos(\theta_2), S2=sin(θ2)S_2 = \sin(\theta_2), C3=cos(ϕ3)C_3 = \cos(\phi_3), and S3=sin(ϕ3)S_3 = \sin(\phi_3). The Euler angles (ϕ1,θ2,ϕ3)(\phi_1, \theta_2, \phi_3) chosen are fairly arbitrarily, as are the signs of rotation matrices given above. This convention was chosen to preserve internal consistency and to produce orientation matrices with sorted directions.

In order to determine the orientation m\tensor{\bm{m}} of a given diffraction pattern, we use a two-step procedure. The first step is to calculate an orientation plan P((ϕ1,θ2),ϕ3,qs)P((\phi_1, \theta_2), \phi_3, q_s) for a given reference crystal. The second step, which is defined in the following section, is to generate a correlogram from each reference crystal, from which we directly determine the correct orientation.

The first two Euler angles ϕ1\phi_1 and θ2\theta_2 represent points on the unit sphere which will become the zone axis of a given orientation. The first step in generating an orientation plan is to select 3 vectors delimiting the extrema of the unique, symmetry-reduced zone axes possible for a given crystal. d shows these boundary vectors for fcc Au, which are given by the directions [001][001], [011][011], and [111][111]. We next choose a sampling rate or angular step size, and generate a grid of zone axes to test. We define a 2D grid of vectors on the unit sphere which span the boundary vectors by using spherical linear interpolation (SLERP) formula defined by . These points with a step size of 22^\circ are shown in d. The rotation matrices which transform the zone axis vector (along the z axis) are given by the matrix inverse of the first two terms in .

We then examine the vector lengths of all non-zero reciprocal lattice points ghkl\bm{g}_{hkl} and find all unique spherical shell radii qsq_s. These radii will become the first dimension of our orientation correlogram, where each radius is assigned one index ss. We loop through all included zone axes, and calculate a polar coordinate representation of the kinematical diffraction patterns.

For each zone axis, the first step to compute the plan is to rotate all structure factor coordinates by the matrix inverse of the first two terms in . Next, we compute the excitation errors sgs_g for all peaks assuming a [0,0,1][0,0,1] projection direction, and the in-plane rotation angle of all peaks γq\gamma_{\bm{q}}. The intensity values of the orientation plan for all qsq_s shells and in-plane rotation values ϕ3\phi_3 are defined using the expression

P0((ϕ1,θ2),ϕ3,qs)={g:g=qs}qsγVgω×max{11δsg2+[mod(ϕ3γg+π,2π)π]2qs2,0},\begin{align*} && P_0((\phi_1, \theta_2), \phi_3, q_s) = \sum_{ \left\{ \bm{g} \,:\, |\bm{g}| = q_s \right\} } {q_s}^\gamma |V_{\bm{g}}|^\omega \times \nonumber \\ && \rm{max} \left\{ 1 - \frac{1}{\delta} \sqrt{ s_{\bm{g}}^2 + \left[ \rm{mod}( \phi_3 - \gamma_{\bm{g}} + \pi, 2 \pi) - \pi \right]^2 {q_s}^2 }, 0 \right\}, \nonumber \end{align*}

where δ\delta is the correlation kernel size, γ\gamma and ω\omega represent the power law scaling for the radial and peak amplitude terms respectively, max(...)\rm{max}(...) is the maximum function, which returns the maximum of its two arguments, mod(...)\rm{mod}(...) is the modulo operator, and the summation includes only those peaks g\bm{g} which belong to a given radial value qsq_s. We have used the combined indexing notation for (ϕ1,θ2)(\phi_1, \theta_2) to indicate that in practice, this dimension of the correlation plan contains all zone axes, and thus the entire array has only 3 dimensions. The correlation kernel size δ\delta defines the azimuthal extent of the correlation signal for each reciprocal lattice point. Note that and are not used for the calculation of orientation plans.

We normalize each zone axis projection using the function

A(ϕ1,θ2)=1ϕ3qsP0((ϕ1,θ2),ϕ3,qs)2,\nonumber A(\phi_1, \theta_2) = \frac{1}{\sqrt{ \sum_{\phi_3} \sum_{q_s} {P_0((\phi_1, \theta_2), \phi_3, q_s)}^2 }},

yielding the final normalized orientation plan

P((ϕ1,θ2),ϕ3,qs)=A(ϕ1,θ2)P0((ϕ1,θ2),ϕ3,qs)P((\phi_1, \theta_2), \phi_3, q_s) = A(\phi_1, \theta_2) P_0((\phi_1, \theta_2), \phi_3, q_s)

By default, we have weighted each term in the orientation plan with the prefactor qsVgq_s |V_g|, i.e.\ setting γ=ω=1\gamma=\omega=1. The qsq_s term gives slightly more weight to higher scattering angles, while the Vg|V_g| term is used to weight the correlation in favour of peaks with higher structure factor amplitudes, which was found to be more reliable than weighting the orientation plan by Vg2|V_g|^2, which weights each peak by its structure factor intensity.

f shows 2D slices of the 3D orientation plan, for the 5 diffraction patterns shown in e. The in-plane rotational symmetry of each radial band is obvious for the low index zone axes, e.g. for the [001][001] orientated crystal, the first row of the corresponding orientation plan consists of four spots which maintains the 4-fold symmetry of the diffraction pattern and can be indexed as [020][020], [200][200], [020][0\overline{2}0] and [200][\overline{2}00]. The final step is to take the 1D Fourier transform along the ϕ3\phi_3 axis in preparation for the Fourier correlation step defined in the next section.

#Correlation Pattern Matching

For each diffraction pattern measurement, we first measure the location and intensity of each Bragg disk by using the template matching procedure outlined by . The result is a set of MM experimental diffraction peaks defined by (qm,γm,Im)(q_m,\gamma_m ,I_m). From these peaks, we calculate the sparse polar diffraction image X(ϕ3,qs)X(\phi_3, q_s) using the expression

X(ϕ3,qs)={qm:qmqs<δ}qmγImω/2×max{11δ(qmqs)2+[mod(ϕ3γm+π,2π)π]2qs2,0}.\begin{align*} && X(\phi_3, q_s) = \sum_{\{ q_m \,: \, |q_m - q_s| < \delta \}} {q_m}^\gamma {I_m}^{\omega/2} \times \rm{max} \left\{ 1 - \right. \\ % \nonumber && \left. \frac{1}{\delta} \sqrt{ (q_m - q_s)^2 + \left[ \rm{mod}( \phi_3 - \gamma_m + \pi, 2 \pi) - \pi \right]^2 {q_s}^2 }, 0 \right\}. \nonumber \end{align*}

By default, we again use prefactors weighted by the peak radius and estimated peak amplitude given by the square root of the measured disk intensities. However, if the dataset being analyzed contains a large number of different sample thicknesses, multiple scattering can cause strong oscillations in the peak amplitude values. As we will see in the simulations below, in these situations the best results may be achieved by setting ω=0\omega=0, i.e. ignoring peak intensity and weighting only by the peak radii. Note that in the diffraction image, the correlation kernel size δ\delta again gives the azimuthal extent of the correlation signal. However, in it also sets the range over which peaks are included in a given radial bin, and the fraction of the intensity assigned to each radial bin.

Next, we calculate the correlation C((ϕ1,θ2),ϕ3)C((\phi_1, \theta_2), \phi_3) of this image with the orientation plan using the expression

C((ϕ1,θ2),ϕ3)=qsF1{F{P((ϕ1,θ2),ϕ3,qs)}F{X(ϕ3,qs)}},\begin{align*} && C((\phi_1, \theta_2), \phi_3) = \sum_{q_s} \ift \left\{ \right. \nonumber \\ && \left. \ft\left\{ {P((\phi_1, \theta_2), \phi_3, q_s)} \right\}^* \ft\left\{ X(\phi_3, q_s) \right\} \right\}, \nonumber \end{align*}

where F\ft and F1\ift are 1D forward and inverse fast Fourier transforms (FFTs) respectively along the ϕ3\phi_3 direction, and the {}^* operator represents taking the complex conjugate. We use this correlation over ϕ3\phi_3 to efficiently calculate the in-plane rotation of the diffraction patterns. The maximum value in the correlogram will ideally correspond to the most probable orientation of the crystal. In order to account for mirror symmetry of the 2D diffraction patterns, we can also compute the correlation

Cmirror((ϕ1,θ2),ϕ3)=qsF1{F{P((ϕ1,θ2),ϕ3,qs)}F{X(ϕ3,qs)}},\begin{align*} && C_{\rm{mirror}}((\phi_1, \theta_2), \phi_3) = \sum_{q_s} \ift \left\{ \right. \nonumber \\ && \left. \ft\left\{ {P((\phi_1, \theta_2), \phi_3, q_s)} \right\}^* \ft\left\{ X(\phi_3, q_s) \right\}^* \right\}, \nonumber \end{align*}

where the mirror operation is accomplished by taking the complex conjugate of F{X(ϕ3,qs)}\ft\left\{X(\phi_3, q_s)\right\}. For each zone axis (ϕ1,θ2)(\phi_1, \theta_2), we take the maximum value of CC and CmirrorC_{\rm{mirror}} in order to account for this symmetry. show 5 output correlograms, for the 5 diffraction patterns shown in e. For each zone axis (ϕ1,θ2)(\phi_1, \theta_2), we have computed the maximum correlation value, which are plotted as a 2D array in g. In each case, the highest value corresponds to the correct orientation.

h shows the correlation values along the ϕ3\phi_3 axis, for the (ϕ1,θ2)(\phi_1, \theta_2) bins with the highest correlation value in g. The symmetry of the correlation values in h reflect the symmetry of the underlying patterns. For the [0,0,1][0,0,1], [0,1,1][0,1,1], and [1,1,1][1,1,1], diffraction patterns, the in-plane angle ϕ3\phi_3 correlation signals have 4-fold, 2-fold, and 6-fold rotational symmetry respectively. By contrast, the asymmetric diffraction patterns with zone axes [1,1,3][1,1,3], and [1,3,5][1,3,5] have only a single best in-plane orientation match.

#Matching of Overlapping Diffraction Patterns

In order to match multiple overlapping crystal signals, we have implemented an iterative detection process. First, we use the above algorithm to determine the best fit orientation for a given pattern. Next, the forward diffraction pattern is calculated for this orientation. We then loop through all experimental peaks, and any within a user-specified deletion radius are removed from the pattern. Peaks which are outside of this radius, but within the correlation kernel size, have their intensities reduced by a factor defined by the linear distance between the experimental and simulated peaks divided by the distance between the correlation kernel size and the deletion radius. Then, the ACOM correlation matching procedure is repeated until the desired number of matches have been found, or no further orientations are found. Note that while we could update the correlation score after peak deletion, we output the original magnitude of the full pattern correlogram in order to accurately calculate the probability of multiple matches.

Examples of alternative orientation plan types in \pyFDSTEM. Fiber texture examples where (a) orientations fully orbit around a single zone axis (the fiber axis), or (b) contain only a symmetry-reduced wedge of zone axes which orbit around a the fiber axis. (c) Examples of orientation plans generated directly from Materials Project entries , using pymatgen symmetries .

#Figure 2:Examples of alternative orientation plan types in py4DSTEM\pyFDSTEM. Fiber texture examples where (a) orientations fully orbit around a single zone axis (the fiber axis), or (b) contain only a symmetry-reduced wedge of zone axes which orbit around a the fiber axis. (c) Examples of orientation plans generated directly from Materials Project entries , using pymatgen symmetries .

#ACOM Integration into py4DSTEM\pyFDSTEM

The ACOM pattern matching described has been implemented into the py4DSTEM\pyFDSTEM python toolkit written by . A typical ACOM workflow starts with using py4DSTEM\pyFDSTEM to import the 4D dataset and one or more images of the vacuum probe. We then use a correlation template matching procedure to find the positions of all diffracted disks at each probe position . We use the correlation intensity of each detected peak as an estimate of the peak's intensity. The resulting set of MM peaks defined by the values (qm,γm,Im)(q_m, \gamma_m, I_m) are stored as a PointList object in py4DSTEM\pyFDSTEM. Because the number of peaks detected at each probe position can vary, we store the full set of all detected peaks in a PointListArray object in py4DSTEM\pyFDSTEM, which provides an interface to the ragged structured numpy data.

Most experimental datasets contain some degree of ellipticity, and the absolute pixel size must be calibrated. We perform these corrections on the set of measured diffraction disks using the py4DSTEM\pyFDSTEM calibration routines defined by . We know that the correlation approach is relatively robust against both ellipticity and small errors in the reciprocal space pixel size. However, precise phase mapping may require us to distinguish between crystals with similar lattice parameters; these experiments will require accurate calibration.

We perform ACOM in py4DSTEM\pyFDSTEM by first creating a Crystal object, either by specifying the atomic basis directly, or by using the pymatgen package to import structural data from crystallographic information files (CIF), or the Materials Project database . The Crystal object is used to calculate the structure factors, and generate an orientation plan. The final step is to use the orientation plan to determine the best match (or matches) for each probe position, from the list of calibrated diffraction peaks. If the sample contains multiple phases, we perform the orientation plan calculation and correlation matching for each unique crystal structure.

In addition to specifying the orientation plan spanning 3 vectors as in , we define additional methods to describe the space of possible orientations. One such example is fiber texture, where we assume the crystals are all orientated near a single zone axis known as the fiber axis, shown in . We can vary the angular range of zone axes included away from the fiber axis as in a, as well as choose the azimuthal range around this axis as in b to account for symmetry around the fiber axis. Alternatively, an "automatic" option is provided, which uses pymatgen to determine the symmetry of the structure and automatically choose the span of symmetrically unique zone axes which should be included in the orientation plan, based on the point group symmetry . This is shown for a selection of different Materials Project database entries in c.

#Simulations of Diffraction Patterns from Thick Samples

One important metric for the performance of an orientation mapping algorithm is how well it performs when the diffraction patterns contain significant amounts of multiple scattering. We have therefore used our ACOM algorithm to measure the orientation of simulated diffraction patterns from samples tilted along many directions, over a wide range of thicknesses. We performed these simulations using the multislice algorithm , and methods defined by and . These methods are implemented in the prismatic simulation code by . The diffraction patterns were generated using a acceleration potential of 300 keV, a 0.5 mrad convergence angle, with real space and reciprocal pixel sizes of 0.05 A˚\rm{\AA} and 0.01 A˚1\rm{\AA}^{-1} respectively, with 4 frozen phonons. In total we have simulated 3750 diffraction patterns from Cu, Ag, and Au fcc crystals, over 25 zone axes ([0,0,1][0,0,1] to [3,4,4][3,4,4] excluding symmetrically redundant reflections) and thicknesses up to 100 nm with a 2 nm step size. These diffraction patterns were generated using the simulation pipeline and database defined by .

#Chemical Synthesis of Twisted AuAgPd Nanowires

The performed synthesis was modified from a known method given by . We prepared the following solutions: 500 mMolar PVP (MW 40,000) in DMF, 50 mMolar HAuCl4 in DMF, 50 mMolar AgNO3 in MilliQ water, and 400 mMolar L-ascorbic acid in MilliQ water. We created the reaction solution in a 4 mL vial (washed 3x with MilliQ water and acetone) by mixing 800 µL DMF, 100 µL PVP, 20 µL HAuCl4, and 20 µL AgNO3. We vortexed the solution, then added 100 µL of L-ascorbic acid solution drop-wise to the mixture while gently swirling. At this point, the color changed from pale yellow to clear. We left the solution at room temperature for 7 days, at which point the solution was light brown/purple. The primary product of this reaction was straight, ultrathin Au-Ag nanowires (2 nm in diameter).

To twist the underlying ultrathin Au-Ag nanowires, we prepared solutions of 875 mMolar L-ascorbic acid and 2 mMolar H2PdCl2 in MilliQ water. In a 4 mL vial (3x washed with MilliQ water/acetone), we added 50 µL of the Au-Ag reacted solution to 640 µL of the L-ascorbic acid solution. Finally, we added 60 µL of the H2PdCl4 solution and allowed the sample to incubate for at least 30 minutes. We purified the reaction solution by centrifuging the product down at 7500 rpm for 4 minutes. We decanted the supernatant, and then rinsed the reaction with MilliQ water 3 times and re-dispersed in MilliQ water. We prepared TEM samples of this material by depositing 10 µL of purified nanowire solution onto 400 mesh formvar/ultrathin carbon grids.

#4D-STEM Experiments with Patterned Apertures

We collected the experimental data using a double aberration-corrected modified FEI Titan 80-300 microscope (the TEAM I instrument at the National Center for Electron Microscopy within Lawrence Berkeley National Laboratory). This microscope is equipped with a Gatan K3 detector and Continuum spectrometer, and was set to collect diffraction patterns integrated over 0.05 seconds, with 4x binning. We used an accelerating voltage of 300 keV, an energy slit of 20 eV, and a spot size of 6. We used a 10 µm bullseye aperture to form the STEM probe in order to improve detection precision of the Bragg disks . We used a convergence semiangle of 2 mrad, with a camera length of 1.05 m. We recorded the experimental dataset using a step size of 5 Å, with a total of 286 and 124 steps in the x and y directions.

Zone axis misorientation as a function of sampling and maximum scattering angle for kinematical simulations. The mean tilt error and number of patterns matched per second are shown inset for each panel.

#Figure 3:Zone axis misorientation as a function of sampling and maximum scattering angle for kinematical simulations. The mean tilt error and number of patterns matched per second are shown inset for each panel.

Dynamical simulated diffraction patterns. (a) Example diffraction patterns for Au oriented to the [011] zone axis for 10-80 nm thick slices. (b) Plots showing the mean zone axis misorientation in degrees as a function of thickness for Cu, Ag, and Au. Each plot shows the errors for correlation prefactors of q_{\rm{s}} |V_{\rm{g}}| (red) and q_{\rm{s}} (blue).

#Figure 4:Dynamical simulated diffraction patterns. (a) Example diffraction patterns for Au oriented to the [011] zone axis for 10-80 nm thick slices. (b) Plots showing the mean zone axis misorientation in degrees as a function of thickness for Cu, Ag, and Au. Each plot shows the errors for correlation prefactors of qsVgq_{\rm{s}} |V_{\rm{g}}| (red) and qsq_{\rm{s}} (blue).

#Results and Discussion

#ACOM of Kinematical Calculated Diffraction Patterns

For the first test of our correlation method, we applied it to the same patterns calculated to generate an orientation plan for fcc Au. Next, we measured the calculation time and angular error between the measured and ground truth zone axes for each pattern. These results are plotted in for 3 different maximum scattering angles kmaxk_{\rm{max}}, and angular sampling of 11^\circ and 22^\circ.

The results in show that the angular error in zone axis orientation is relatively insensitive to the angular sampling. However, the angular error drops by a factor of 10 from approximately 3\approx 3^\circ to 0.3\approx 0.3^\circ when increasing the maximum scattering angle included from kmax=1A˚1k_{\rm{max}} = 1 \, {\rm{\AA}}^{-1} to 1.5A˚11.5 \, {\rm{\AA}}^{-1}, and by another factor of 2-3 when increasing kmaxk_{\rm{max}} to 2A˚12 \, {\rm{\AA}}^{-1}. This is unsurprising, as examining e shows that there is a large increase in the number of visible Bragg spots outside of kmax=1A˚1k_{\rm{max}} = 1 \, {\rm{\AA}}^{-1}, and because Bragg disks at higher scattering angles provide better angular precision relative to low kk disks. This result emphasizes the importance of collecting as wide of angular range as possible when performing orientation matching of 4D-STEM data.

The inset calculation times reported are for the single-threaded ACOM implementation in py4DSTEM\pyFDSTEM, running in Anaconda on a laptop with an Intel Core i7-10875H processor, running at 2.30 GHz. The calculation times can be increased by an order of magnitude or more when running in parallel, or by using a GPU to perform the matrix multiplication and Fourier transform steps.

#ACOM of Dynamical Simulated Diffraction Patterns

In diffraction experiments using thick specimens, the electron beam can scatter multiple times, a phenomenon known as dynamical diffraction. This effect is especially pronounced in diffraction experiments along low index zone axes, where the diffracted peak intensities oscillate a function of thickness. In order to test the effect of oscillating peak intensities on our ACOM method, we have simulated diffraction patterns for Cu, Ag, and Au fcc crystals, along multiple zone axes. Some example diffraction patterns for the [011] zone axis of Au are plotted in a. We see that all diffraction spots have intensities which oscillate multiple times as a function of thickness.

We performed ACOM by generating orientation plans with an angular sampling of 22^\circ, a correlation kernel size of 0.08 A˚1\rm{\AA}^{-1}, and maximum scattering angles of kmaxk_{\rm{max}} = 1.0, 1.5, and 2.0 A˚1\rm{\AA}^{-1}. We kept the radial prefactor of weighting set to γ\gamma = 1, and tested peak amplitude prefactors of ω\omega = 1.0, 0.5 and 0.0. The average zone axis angular misorientation as a function of thickness is plotted in b. In total, we performed orientation matching on 3750 diffraction patterns, and a total of 33750 correlation matches on a workstation with an AMD Ryzen Threadripper 3960X CPU (2.2 GHz, baseclock). The typical number of patterns matched per second were of ~80-90, 45-55 and 25-30 patterns/s for kmaxk_{\rm{max}} values of 1.0, 1.5, and 2.0 A˚1\rm{\AA}^{-1} respectively.

As expected, the errors are higher than those achieved under kinematic conditions, and the trend for smaller errors with larger kmaxk_{\rm{max}} is also preserved (mean errors of 7.25^\circ, 3.09^\circ and 1.39^\circ for kmaxk_{\rm{max}} values 1.0, 1.5 and 2.0 A˚1\rm{\AA}^{-1} respectively, γ\gamma = 1, ω\omega = 0.25). We did not observe any dependence of the orientation accuracy on the simulation thickness. Despite the correlation prefactor Vg|Vg| performing well for the examples shown in , for the dynamical diffraction simulations along zone axes it was out-performed by prefactors of both Vg\sqrt{|V_g|} (ω\omega = 0.5) and omitting the peak amplitude prefactor altogether (ω\omega = 0). We therefore suggest that when mapping samples with a large range of thicknesses, or many crystals aligned to low index zone axes, the position of the diffracted peaks is significantly more important than their amplitudes or intensities. One possible method to increase the accuracy while using higher amplitude prefactors is to perform an experiment which recovers more kinematical values for the diffracted peak intensities, for example by precessing the electron beam when recording diffraction patterns . We note that there is likely no global optimal choice of orientation mapping hyperparameters for all materials and thicknesses, and this may be a worthwhile topic for future investigations.

#4D-STEM ACOM of Twisted AuAgPd Nanowires

We have tested our ACOM algorithm with a 4D-STEM dataset collected for an AuAgPd nanowires. An image of the vacuum bullseye STEM probe is shown in a. For each detector pixel, we have calculated the maximum value across all STEM probe positions to generate a maximum diffraction pattern, shown in b. The beamstop used to block the center beam is visible, as well as various crystalline diffraction rings out to approximately 1.4 A˚1{\rm{\AA}}^{-1}.

4D-STEM scan of twisted polycrystalline AuAgPd nanowires. (a) Diffraction image of probe over vacuum, showing bullseye pattern. (b) Maximum of each pixel in diffraction space over all probe positions. (c) Histogram of all peak locations detected by correlation in \pyFDSTEM of (a) with each pattern included in (b). (d) HAADF-STEM image of the sample. (e) 1D histogram of scattering vectors, with fcc AuAg inverse plane spacings overlaid.

#Figure 5:4D-STEM scan of twisted polycrystalline AuAgPd nanowires. (a) Diffraction image of probe over vacuum, showing bullseye pattern. (b) Maximum of each pixel in diffraction space over all probe positions. (c) Histogram of all peak locations detected by correlation in py4DSTEM\pyFDSTEM of (a) with each pattern included in (b). (d) HAADF-STEM image of the sample. (e) 1D histogram of scattering vectors, with fcc AuAg inverse plane spacings overlaid.

After performing the correlation peak finding algorithm in py4DSTEM\pyFDSTEM, we have an estimated position and intensity of all detected Bragg peaks. A 2D histogram of these peaks, known as a Bragg vector map, is plotted in c. Sharp polycrystalline diffraction rings are clearly visible, as well as false positives generated by the beamstop edge. These false positives were manually removed by using a mask generated from an image of the beamstop. A high angle annular dark field (HAADF) image was simultaneously recorded during the 4D-STEM data collection, which is shown in d.

The final experimental pre-processing steps are to calibrate the diffraction pattern center, the elliptical distortions, and the absolute pixel size. We performed these steps by fitting an ellipse to the (022)(022) diffraction ring, and by assuming a lattice constant of 4.08A˚4.08 \, \rm{\AA}, corresponding to the fcc Au structure . This process is explained in more detail by . We assumed that the Ag lattice constant is similar to that of Au. Despite the presence of Pd in the nanowires, there was no significant presence of secondary grains corresponding to the smaller lattice of fcc Pd grains. An intensity histogram of the corrected Bragg peak scattering angles are shown in e. We have overlaid the 5 smallest scattering angles of Au on e to show the accuracy of the correction.

Orientation mapping of polycrystalline AuAgPd nanowires. (a) Total of measured correlation signal for each probe position. (b) Estimated number of patterns indexed for each probe position. (c) Example of 2 orientations indexed from a single diffraction pattern, collected at the position indicated by the arrow shown in (b). (d) Orientation maps of the 3 highest correlation signals for each probe position. Legend shown above.

#Figure 6:Orientation mapping of polycrystalline AuAgPd nanowires. (a) Total of measured correlation signal for each probe position. (b) Estimated number of patterns indexed for each probe position. (c) Example of 2 orientations indexed from a single diffraction pattern, collected at the position indicated by the arrow shown in (b). (d) Orientation maps of the 3 highest correlation signals for each probe position. Legend shown above.

We have performed ACOM on the AuAgPd nanowire sample, with the results shown in shown for up to 3 matches for each diffraction pattern. For each probe position, the sum of the maximum detected correlation signals for up to three matches are shown in a. The structure is in good agreement with , though with additional modulations due to some grains generating more diffraction signal than others. Using a correlation intensity threshold of 0.5, we have plotted the number of matching patterns in b. The threshold of 0.5 was arbitrary chosen as a lower bound for a potential match. Examples of 2 matches to a single diffraction pattern are plotted in c. In this figure the correlation score for the first matched pattern was higher than the second. The second match found shows some deformation between the measured and simulated Bragg peak positions, and matches less peaks. It therefore produces a lower correlation score, which can be used to threshold the results as in d.

d shows the 3D orientations for all probe positions, with the 3 best matches shown. Each image is masked by the total correlation signal, so that low correlation values are colored black. Almost every diffraction pattern with Bragg disks detected was indexed for at least one orientation with high confidence. Additionally, the patterns are very consistent, with a large number of adjacent probe positions recording the same orientation. Some secondary grains are also clearly visible in the second-best match, while very few patterns have been assigned a third match with high confidence.

Orientation analysis of grains in AuAgPd nanowires. (a) Crystal grains, with in-plane (111) planes colored by orientation. (b) (111) planes shared by two overlapping grains.

#Figure 7:Orientation analysis of grains in AuAgPd nanowires. (a) Crystal grains, with in-plane (111) planes colored by orientation. (b) (111) planes shared by two overlapping grains.

In order to investigate the grain organization of the AuAgPd nanowires, we have performed clustering analysis on the orientation maps. Grains with similar orientations have been clustered together by looping through each probe position and comparing its orientation to its neighbors. Grains with at least 10 contiguous probe positions are shown in a. (111)(111) planes which lie in the image plane are overlaid onto the grain strucure, colored by their orientation. Confirming our observations in d, only a few grains with substantial overlap were reliably identified. This might be due to the low thickness of the sample (only a single grain along the beam direction), some grains not being oriented close enough to a zone axis to be detected, or multiple scattering deviations in the diffracted signal. There is a noticable bias in the orientation of the (111)(111) planes, which tend to be oriented horizontally near the growth direction of the nanowires.

One hypothesis for the growth mode of these twisted nanowires is that adjacent grains are connected by (111)(111) twin planes, forming local helical structures to give the observed twisted structures. To test this hypothesis, we determined the position of (111)(111) planes from a which are shared by two overlapping grains. b shows the location of these shared (111)(111) planes (with plane normal differences below 88^\circ), colored by the normal vector of the plane. Many shared (111)(111) planes were detected, most with normal vectors aligned to the wire growth direction. These observations support the hypothesis that these nanowires are composed of grains connected by (111)(111) twin planes.

These experimental observations demonstrate the efficacy of our ACOM method. In order to improve these results, we will need to collect diffraction data with a wider angular range. This can be achieved by using precession electron diffraction , multibeam electron diffraction , or by tilting the sample or beam and recording multiple 4D-STEM datasets .


We have introduced an efficient and accurate method to perform automated crystal orientation mapping, using a sparse correlation matching procedure. We have implemented our methods into the open source py4DSTEM\pyFDSTEM toolkit, and demonstrated the accuracy of our method using simulated diffraction patterns. We also applied ACOM to an experimental scan of a complex helical polycrystalline nanowire, where we were able to identify shared twin planes between adjacent grains which may be responsible for the twisted helical geometry. All of our methods have been made freely available to the microscopy community as open source codes. We believe that our implementation of ACOM is efficient and accurate enough to be Incorporated into automated online TEM software . In the future, we will improve our ACOM method using machine learning methods , and we will extend our ACOM methods to include multibeam electron diffraction experiments .

#Source Code and Data Availability

All code used in this manuscript is available on the py4DSTEM GitHub repository, and the tutorial notebooks are available on the py4DSTEM tutorial repository. All simulated and experimental 4D-STEM datasets are available at [links will be added after publication].


We thank Karen Bustillo for helpful discussions. CO acknowledges support of a US Department of Energy Early Career Research Award. SEZ was supported by the National Science Foundation under STROBE Grant no. DMR 1548924. AB, BHS, and py4DSTEM\pyFDSTEM development are supported by the Toyota Research Institute. AR is supported by the 4D Data Distillery project, funded by the US Department of Energy. Work at the Molecular Foundry was supported by the Office of Science, Office of Basic Energy Sciences, of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. This research used resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility located at Lawrence Berkeley National Laboratory, operated under Contract No. DE-AC02-05CH11231.

  1. Borchardt-Ott, W. (2011). Crystallography: an introduction. Springer Science & Business Media.
  2. Dederichs, P., Lehmann, C., Schober, H., Scholz, A., & Zeller, R. (1978). Lattice theory of point defects. Journal of Nuclear Materials, 69, 176–199.
  3. LeSar, R. (2014). Simulations of dislocation structure and response. Annu. Rev. Condens. Matter Phys., 5(1), 375–407.
  4. Tang, M., Carter, W. C., & Cannon, R. M. (2006). Diffuse interface model for structural transitions of grain boundaries. Physical Review B, 73(2), 024102.
  5. Janssen, G. (2007). Stress and strain in polycrystalline thin films. Thin Solid Films, 515(17), 6654–6664.