## research papers

## of high-energy electrons by light-atom structures: a multiple interpretation

^{a}STFC, Rutherford Appleton Laboratory, Didcot, OX11 0FA, United Kingdom, and ^{b}CCP4, Research Complex at Harwell, Rutherford Appleton Laboratory, Didcot, OX11 0FA, United Kingdom^{*}Correspondence e-mail: eugene.krissinel@stfc.ac.uk

Because of the strong electron–atom interaction, the kinematic theory of diffraction cannot be used to describe the scattering of electrons by an assembly of atoms due to the strong

that needs to be taken into account. In this paper, the scattering of high-energy electrons by a regular array of light atoms is solved exactly by applying the T-matrix formalism to the corresponding Schrödinger's equation in spherical coordinates. The independent atom model is used, where each atom is represented by a sphere with an effective constant potential. The validity of the approximation and the phase grating approximation, assumed by the popular multislice method, is discussed, and an alternative interpretation of multiple scattering is proposed and compared with existing interpretations.Keywords: high-energy electron diffraction; T-matrix; multiple scattering; independent atom approximation.

### 1. Introduction

#### 1.1. Motivation

The theory of ; Cowley, 1995), very shortly after the first experimental demonstration of electron diffraction. Over the 20th century, several theories of electron diffraction have been extensively applied in various areas of solid-state physics. In particular, the multiple scattering theory (MST) (Korringa, 1947, 1994; Kohn & Rostoker, 1954; Dederichs, 1971) has proven to be very efficient for describing electronic properties of matter. Besides, MST provides an intuitive picture of Partial wave scattering theory (Schiff, 1955) provides another description of multiple scattering with the possibility to take into account (Howie, 1963) and loss of coherency (Howie, 2014). Other less commonly known attempts at an intuitive depiction of have been proposed to explain electron channelling in crystals, for example by VanDyck & Op de Beeck (1996). must be taken into account when considering electron–atom interactions, since even at the very high electron energies commonly used in modern transmission electron microscopes, the interaction is so strong that, theoretically, the kinematic approximation is not valid for ideal crystals thicker than 10–15 nm (Hirsch *et al.*, 1965; Glaeser & Downing, 1993; Subramanian *et al.*, 2015). In practice, crystal growth cannot be controlled to such a degree of accuracy and, usually, nanocrystals of organic compounds have a size on the order of hundreds of nanometres. This is a challenging aspect of high-energy electron diffraction in crystallography, as it significantly complicates the process. However, structures have been successfully determined from ED (electron diffraction) data using the standard kinematic theory (Gemmi *et al.*, 2019; Nannenga *et al.*, 2014; Nannenga & Gonen, 2019). Although dynamical (Palatinus *et al.*, 2013) usually leads to better intensity predictions (Gemmi *et al.*, 2019), the agreement between theory and experiment is still significantly worse than that obtained for X-ray data (Oleynikov *et al.*, 2007; Klar *et al.*, 2021).

Although there is no quantitative prediction of the effects of

on the ability to solve a it may be speculated that may be attenuated by other effects such as a lack of coherency caused by solvent scattering or crystal defects. Therefore, more accurate theoretical models are needed for getting better results from electron diffraction experiments.#### 1.2. State of the art

The theory of multiple scattering has a long and rich history, including modern theoretical developments and implementations (Sébilleau *et al.*, 2017; Zabloudil *et al.*, 2005). However, in the context of high-energy electron diffraction, the multislice (MS) (Cowley & Moodie, 1957) approach and Bloch-wave (BW) (Bethe, 1928; Metherell & Fisher, 1969) approach are the most popular methods for simulating diffraction patterns by crystals.

MS is particularly well suited for solving large problems as it involves successive convolutions, which can be very efficiently computed with the fast Fourier transform (FFT) (Ishizuka & Uyeda, 1977). However, in order to avoid aliasing, transverse periodic boundary conditions must be met. This is not always possible if the crystal is in an arbitrary orientation as is the case in typical continuous-rotation electron diffraction experiments. A small beam tilt can also be used as an option to remedy this aspect but only extends reliably to a maximum 3–6° (Ishizuka, 1982; Chen *et al.*, 1997). Simulations with arbitrary orientations could also be performed by considering all unit cells in a crystal, while adding zero padding at its ends. Even though modelling a full crystal can quickly become intractable computationally, the use of a smooth envelope function may improve this commonly admitted limitation of MS (Kirkland, 2019). A few very efficient MS implementations have been reported to date (Ophus, 2017). Some of them are targeted at convergent-beam electron diffraction (CBED), while others have special features such as inclusion of (Allen *et al.*, 2015). Yet, an approach combining full modelling capabilities with an acceptable computational efficiency, designed specifically for the continuous-rotation electron diffraction of large organic structures, does not seem to be available. Such an approach would indeed be of significant value for macromolecular from ED data.

On the other hand, the BW method can simulate ED for structures in arbitrary orientations but does not scale well with the structure size. Although a few rather efficient BW implementations (Zuo & Weickenmeier, 1995) have been proposed to include a large number of beams, BW can hardly be applied to very large structures since the becomes prohibitively large. Due to the unfavourable *O*(*N*^{3}) complexity of matrix diagonalization, where *N* is the number of beams, the method is hardly applicable in practice to the determination of macromolecular structures where it would have to be done numerous times. Non-periodic structures, defects and solvent scattering are also hard to model with this method.

In other areas of physics, specialized multiple scattering approaches have also been developed. In particular, the T-matrix approach has been extensively applied to wave propagation related problems such as electromagnetics (Hamid *et al.*, 1990*a*,*b*; Eremin *et al.*, 1995), optics (Moine & Stout, 2005) and acoustics (Silva *et al.*, 2012; Godin, 2011). Although it does not compete with MS from a computational standpoint, this approach benefits from providing an exact solution to Schrödinger's equation for an ensemble of spherically symmetric atomic potentials in the independent atom model (IAM) approximation. Besides, *ab* *initio* real-space multiple scattering calculations have also been developed in X-ray photon emission spectroscopy with elaborate potentials including many-body effects (Rehr *et al.*, 2009). Similar approaches based on the T-matrix also exist in electron energy-loss spectroscopy (Sébilleau *et al.*, 2006).

#### 1.3. Contribution and outline

The purpose of this paper is to adapt the T-matrix formalism specifically to the case of scattering of fast electrons by light-atom structures. In particular, we present an intuitive picture of multiple scattering in the

approximation. We propose a comparison with the MS interpretation of multiple scattering. In addition, we discuss the validity of the approximation and the phase grating approximation used in the MS approach. Finally, we draw conclusions and outline extensions of the proposed approach in order to account for incoherent inelastic electron scattering.### 2. Theory

The scattering of fast electrons by an atomic structure is described by the wavefunction Ψ, obeying the following Schrödinger's equation (Kirkland, 2019):

where is Planck's constant, the mass of the electron, *e* the the spatially varying electrostatic potential created by the atoms, and *E* the electron energy.

In an ED experiment, , which corresponds to the continuum, non-quantized state of the system (Chuang, 1996). Therefore, we solve (1) with *E* as a parameter, equal to the energy of incident electrons.

The total wavefunction Ψ is described as the sum of an incident wave and a scattered wave :

The incident wave is assumed to be known, *i.e.* a plane wave or a Gaussian wave for example. The T-matrix aims at determining the scattered wave by the sample.

#### 2.1. T-matrix formulation

In its standard form, the T-matrix approach solves for the case where the incident wave is described by a plane wave of wavenumber (optics convention), and the electrostatic potential is modelled by a uniform constant inside non-overlapping spheres. The constant potential and radius of the spheres depend on the atom type.

Although this is not an accurate representation of the actual potential profile *V*(*r*), using an appropriate amplitude for the constant potential *V*_{0} should provide similar values of scattering cross sections. The scattering being a reliable figure of merit of the strength of an interaction, one can assume that a suitable choice of radius and constant potential should therefore provide a faithful depiction of the extent of multiple scattering in the structures of interest.

The setup is shown in Fig. 1. This formulation is well established (Hamid *et al.*, 1990*a*,*b*; Eremin *et al.*, 1995) and the theory is outlined here for the purpose of introducing the forward multiple scattering approximation. In the T-matrix approach, the electrostatic potential is assumed constant inside different regions of the 3D space. In each region, the problem is reduced to the Helmholtz equation in spherical coordinates:

where is the constant positive potential inside sphere *p* of radius *a*_{p} centred at , *k*_{p} the wavenumber inside the sphere and *n*_{p} is sometimes referred to as the by analogy with optics. In the remainder, *n*_{p} will be referred to as the potential strength for clarity purposes. The scattered wavefunction can be decomposed into the scattered wavefunctions of individual atom spheres,

where *N* is the number of spheres. The *p*-sphere scattered wave can itself be partitioned into a part inside sphere *p* noted and a part outside sphere *p* noted . The expressions for and are detailed as

where , *j*_{l} and *h*_{l}^{(1)} are the spherical Bessel and Hankel functions of the first kind, respectively, *Y*_{l}^{m} are the spherical harmonics of order *l* and azimuthal order *m*. Note that these equations are expressed in the reference frame of each sphere *p*, hence the use of variable .

The unknown coefficients , are found by imposing the continuity of the wavefunction and its radial derivative at the surface of each sphere *p*. After some mathematical manipulations (see Appendix *A*), this results in the following linear system of equations:

where **I** is the identity matrix, **A** the unknown vector of coefficients, **T** the cross-coupling matrix and **L** the known incident wave which appears on the right-hand side of (6).

Note that, in this formulation, the unknown coefficients, found from the continuity interface conditions, will model the response of each individual sphere to both the incident wave and the waves scattered by the other atoms. This model therefore does account for multiple scattering.

#### 2.2. Far field and scattering cross section

In electron crystallography, diffraction patterns are recorded that correspond to the intensity of the scattering amplitude profile for various crystal orientations. A diffraction pattern is obtained in the far-field diffraction regime, which arises when the shape of the angular radiation scattering amplitude no longer depends on the distance to the scattering specimen. Using the asymptotic behaviour *h*_{l}^{(1)}(*k*_{0}*r*_{p}) and since , , the far-field scattering amplitude from each individual sphere *p* can be written as

The total far-field scattering amplitude is the sum of the scattered field contributions from all individual spheres. Note that the scattered field from each individual sphere already accounts for the scattering between the spheres through the coefficients *a*_{p} and *b*_{p}. Since, in the far field, ,

The normalized differential scattering

can be obtained aswhere we have used .

### 3. Real-space forward multiple scattering picture

#### 3.1. T-matrix approximations

Equation (6) is a convenient way to represent the system as it readily identifies **L** as the solution to the uncoupled (**T** = 0) problem. If **T** = 0, the system can be broken down separately into the scattering problems for each individual sphere, for which the solutions are immediately identified as , . These are indeed the well known analytical solutions of by a soft sphere (Balanis, 1989).

The cross-coupling matrix **T** accounts for multiple scattering effects. If **A** is written as , then **T** is a matrix with only off-diagonal components:

where **T**_{pq} represents the scattering from sphere *p* due to the scattering from sphere *q*. By considering all components , , we assume that the scattering from sphere *p* affects scattering from sphere *q* and vice versa. Matrix (10) is completely filled with nonzero entries apart from its diagonal components. It is therefore not a sparse matrix, which implies quadratically growing memory requirements. Moreover, since full intersphere scattering is considered, the solution of equation (6) requires matrix inversion which is computationally expensive.

However, as detailed below in Section 4.2, backscattering can be neglected for very fast electrons, which is known as the approximation. This results in **T** being lower triangular if the spheres are sorted in ascending order along . In the case of atoms lying in the same coordinate plane, *i.e.* slice, the atoms can be sorted according to their transverse coordinates. Then, as 90° scattering is neglected in the approximation, the T-matrix remains lower triangular. Inversion of the lower triangular matrix is computationally tractable and calculations can be performed sequentially one slice after the other. Therefore, the approximation converts the implicit self-consistent T-matrix scheme into a explicit scheme similar to the approach taken in the MS method.

#### 3.2. Multiple scattering approximations

Since represents single scattering, we can establish that accounts for secondary scattering. Similarly, outward scattering amplitudes from sphere *p* can be written as

where the first term accounts for kinematic scattering, the second term for secondary scattering, the third term for tertiary scattering, and so on. This is a development similar to the Korringa–Kohn–Rostoker (KKR) theory of multiple scattering (Korringa, 1947, 1994; Kohn & Rostoker, 1954). This forward multiple scattering picture is illustrated in Fig. 2 for the case of a two-atom system. From a computational point of view, equation (11) may offer an advantage over the full approximation if only a few *n*-time scattering terms are necessary. Indeed, since each term in the expansion depends on precomputed coefficients , computation of the coefficients for the *n*th scattering term can be massively parallelized for each atom.

#### 3.3. Scattering probabilities

The probability of an electron being scattered elastically *n* times can be estimated classically from ballistic arguments using a continuous model of matter. It can be established (Egerton, 2011) that the probability of an electron being elastically scattered *n* times after passing through an amorphous sample of thickness *z* follows a , where is the elastic It is paramount to note that this classical expression models not including interference effects between multiple scattering events. This model is therefore valid for classical particles not exhibiting particle–wave duality behaviour. In quantum systems, it is also possible to consider if is strongly present. Indeed, not only induces a transfer of energy between scatterers but also a loss of coherency through phase shifts (Howie, 2014). As a result, quantum systems, where is predominant, would scatter incoherently. For example, may be due to electron–electron interactions, solvent or impurity scattering, as well as thermal diffuse scattering.

In this section, we suggest a way to consider probabilities of multiple coherent

events. It will be shown that while multiple scattering amplitudes can be considered, it may not be possible to identify corresponding probabilities without any ambiguity precisely because of interference.Using a wave-like description, the probability of an electron undergoing a scattering event can be determined from the scattering *S* be the area over which the specimen is illuminated.

Let *z* be the thickness of the specimen. The incident electron plane wave can be normalized by so that it integrates over the interaction volume *V* = *Sz* to a single electron. Then, the flow of electrons per unit time per unit area (Vainshtein, 1964) is defined as , where *v*_{0} is the speed of the incident electrons.

By definition of the *J*_{0} and ,

where *S* can for example be the area of the aperture in (SAED). When performing a MS simulation, *S* is the area of the simulated domain (transverse if performing a periodic simulation).

Here, we define *f*_{n} as the complex scattering amplitudes obtained by considering scattering from the wave being scattered *n* times. It can readily be established that

We can then define as the probability of an electron being scattered *n* times. The choice of *n*-time scattering is arbitrary precisely because of the cross-coupling interference terms *f*_{n}*f*_{m}^{*}. It has been chosen in such a way that when more *n*-scattering terms are taken into account, probabilities of higher *n*-level scattering events are transferred from the lower *n*-level scattering events. In the case of an application to two-beam scattering theory, the intensity of the forward and scattered beams oscillates as the thickness of the sample grows. The scattered beam is zero when the forward beam has been scattered an even number of times and maximum when the forward beam has been scattered an odd number of times (Cowley & Moodie, 1957). Our definition of scattering probabilities should reflect this behaviour.

The probability of an electron not being scattered is naturally , where . Note that, since , the probability of scattering is necessarily less than 1.

For a very large number of scatterers, *N*, regularly spaced by distance *dz*, we can define the average . The corresponding scattering probability is then *P*(*dz*) = where is the density (since *Sdz* contains only one scatterer) and is the This is consistent with the definition used by Latychevskaia & Abrahams (2019).

#### 3.4. First Born approximation and kinematic scattering

Let us assume normal incidence with , so that *e*_{p} = exp(*jk*_{0}*d*_{p}), and note that the can be represented by **q** = in (or transfer vector momentum space), with being the wavenumber in crystallographic convention. Keeping only the kinematic term in equation (11) and inserting in the far field in equation (8) give

where *f*_{p}^{(e)} is the scattering amplitude of sphere *p* in the first Born approximation, also known as the form factor. Equation (13) corresponds to the standard kinematic formula for the traditionally used in crystallography. Note that the kinematic approximation is nothing else than the first Born approximation applied to the whole assembly of atoms.

It will be shown in Section 4.2 that the first Born approximation holds in the case of single scatterers, but fails in application to multi-atom systems. This means that the kinematic approximation is not generally valid.

#### 3.5. Multiple scattering in the MS approach

In MS, the

approximation is used and the potential discretized in slices. The wavefunction is propagated from one slice to the other using a Fresnel propagator.A multiple scattering interpretation of the MS approach has been proposed (Cowley & Moodie, 1957), which, although analogous to the one presented above, differs in that it is stated in More details on the differences between the two interpretations can be found in Appendix *B*.

### 4. Application and results

Very efficient open-source packages exist for T-matrix calculations, but they are targeted at applications in electromagnetics (Egel *et al.*, 2017; Kottke, 2020). There are also packages targeted at photon and electron spectroscopy and related to the T-matrix, such as *msSpec* (Sébilleau *et al.*, 2011) and *FEFF9* (Rehr *et al.*, 2010). Here, an open-source package, adapted to high-energy electron diffraction and including the approximation, has been developed and made available (Drevon, 2021).

#### 4.1. Validity of the implementation

In practice, the size of the matrix in equation (6) has to be truncated to a maximum integer order due to computer memory limitations.

A rule for obtaining accurate results is to use an integer larger than the maximum value of the normalized radius *ka*, which we call *ka*_{max} = max_{p} *ka*_{p} where *a*_{p} is the radius of the spheres and *k* the wavenumber.

Besides, the translational addition theorem (Dufva *et al.*, 2008) used for expressing the scattered field of any sphere in the reference frame of another sphere introduces an approximation of the translated spherical Hankel functions. This theorem is also more commonly known as the structure constant expansion or Kasterin expansion. The computational accuracy of this translation decreases from the centre at which this expansion is written. This is analogous to a Taylor expansion, for which the accuracy away from the point at which it is expressed increases with increasing expansion order. Note that the distances between the spheres do not affect the accuracy of this translation operation and should not affect the choice of . Once again, the radius of the spheres determines the choice of . This is illustrated in Fig. 3(*a*), where the error between *h*_{l}^{(1)} *Y*_{42}, computed at the origin, using the translational addition theorem with and , is displayed in log scale.

The choice of is best evaluated by assessing the continuity of the wavefunction at the surface of the spheres as shown in Fig. 3(*b*) for *N* = 4 and *ka*_{max} = 4. This is also used to validate the correctness of the implementation since it can be seen that machine accuracy can be reached when increasing the order of the expansion.

#### 4.2. Validity of and phase grating approximations for light atoms

In the case of very fast electrons typically used in transmission electron microscopes, *E* = 50–300 keV. Inclusion of relativistic effects results in a wavelength λ = 0.02508 Å at 200 keV, which will be assumed from now on unless stated otherwise.

In the IAM approximation, the Coulomb potential is created by the charge of the nucleus and its electron cloud. It is typically fitted with a sum of three screened Coulomb potentials and three Gaussian terms (Kirkland, 2019) which results in the potential shown in Fig. 4(*a*).

The solution to Schrödinger's equation in such a potential can only be solved perturbatively (Müller, 1965) and is beyond the scope of this study. However, indicative values for *k*_{p} and *ka*_{p} can be used with a multi-shell representation as shown in the form of blue patches in Fig. 4(*a*). Although the range of the screened Coulomb potential is theoretically infinite, it can be truncated to radius *ka* for most practical purposes. This Coulomb potential truncation is commonly known as the muffin-tin model. The number of concentric spheres to use has mainly an effect on the large-angle representation of the form factor which does not impose severe constraints in a low-angle approximation. In Fig. 4(*b*), the multi-shell scattering amplitudes calculated in the Born approximation are shown for increasing values of truncation radius. The figure only shows a satisfactory agreement with the electron diffraction scattering factors for normalized radius as large as *ka* = 400. Even though the potential *V*(*r*) is extremely small at such a large value of the radius *ka* = 400, those shells are required to account for the proper low-angle representation of the form factor.

As mentioned above, ED simulations based on the T-matrix approach are very expensive computationally at *ka* = 400, because higher-order terms of series in equations (4), (5) would need to be included. However, the medium truncation range should be sufficient for providing a reasonable picture of dynamical scattering. Fig. 5(*a*) shows the total scattering of a single sphere with increasing radius for a range of values of potential strength *n*_{p}. The dark blue curve shows the locations of the spherical shell for carbon. It becomes almost flat for the parameter set ( *ka* = 30, *n*_{p} = 1.001), from which it keeps increasing very slightly to reach the asymptotic value . This value is almost identical to the average elastic for real carbon in the Born approximation (Latychevskaia & Abrahams, 2019) which confirms that the chosen parameters capture the strength of incident electron–carbon interaction. Therefore, even though the chosen parameter set is a coarse representation of the potential profile of a real carbon atom, it should provide a reasonably faithful estimate of the extent of multiple scattering of real carbon-based samples, which is the main objective of this study.

Fig. 5(*b*) shows the shape of the far-field amplitudes for two sets of *ka*, *n*_{p} values. The first one is for the extreme value of *ka* = 3 and *n*_{p} = 1.1 while the second corresponds to one of the carbon atoms. While differences are obvious for the first case, the Born approximation is very accurate for the carbon atoms, although some minor differences appear at low angles where the phase grating approximation provides some improvements too.

Fig. 6(*a*) shows a (*ka*,*n*_{p}) map of the difference between the exact coefficients of the two-body problem and those calculated using the kinematic approximation. This error is calculated as

Fig. 6(*b*) shows the same thing as a result of using the forward approximations.

For this case, a value of *kd* = 3*ka* is used, which is close to interatomic distances. Overall, it is clear that the approximation is very accurate over all ranges of parameter sets for carbon atoms, which is expected since there is very little backscattering beyond 90°. On the other hand, the kinematic approximation does not appear quite as good, even for this, a mere two-atom problem. Although not shown here, both the kinematic and approximations tend to work slightly better with increasing distances *kd*, since coupling between the spheres gets reduced and, therefore, is less likely to affect scattering from the other spheres. Using low values of *n*_{p} results in an overall good approximation of both the uncoupled and approximation. This is an anticipated result, since for weak potentials the kinematic approximation works better. The uncoupled approximation improves with decreasing radii, since a small *ka* results in a small scattering On the other hand, the approximation improves with larger values of *ka*, which make backward scattering less likely.

Fig. 7(*a*) shows the scattering coefficients in complex space for the given two-body problem, where *n*_{p} = 1.01, *kd* = 3*ka*, while increasing radius *ka* (up to *ka* = 15) for a selected set of orders *l*. Only the coefficients for sphere *p* = 2 are shown since it is the sphere most affected by the approximations. Here, only *m* = 0 coefficients are used, because in the case of planar illumination, this configuration has azimuthal symmetry. The coefficients are computed using the full T-matrix, the approximation and the kinematic approximation. In order to get this picture, the coefficients were divided by the phase factor at sphere *p*, *i.e.* . The arms of the spiral are rotated by from one order to the other as a result of the *j*^{l} factor in the spherical expansion of the incident plane wave. All the coefficients increase with increasing *ka* as a result of the increasing scattering with increasing *ka*. For a given *ka*, the strongest contributing term orders are around *l* = *ka* as a direct manifestation of the convergence behaviour presented above. It can be observed that the kinematic approximation error increases with the radius while the forward approximation is very accurate across all radii.

Fig. 7(*b*) shows a comparison of the computed far-field scattering amplitudes for *N* = 2, *ka* = 30, *n*_{p} = 1.01, *kd* = 2*ka* using the full T-matrix, the forward approximation and the kinematic approximation. The peaks and valleys are both due to the single scattering profile of the constant sphere and the interferometric path length exp(*jkd*). It is possible to see that the coupled problem has a slight averaging effect over the kinematic pattern, a well known feature of dynamical diffraction.

#### 4.3. Successive multiple scattering approximations

Finally, we consider an array of *N* identical scatterers regularly spaced by *kd* = 3*ka* which correlates with lengths in molecules.

Here, we consider the successive multiple scattering approach under normal illumination. This is an extreme case where strong dynamical scattering is expected.

Figs. 8(*a*) and 8(*b*) show the evolution of the error [equation (14)] on scattering amplitude coefficients for an increasing number of spheres while also varying the number of successive approximations while keeping parameters *ka* = 7, *kdka* = 3 fixed.

The ). Overall, the successive approximation schemes naturally converge to the forward approximation for sufficiently large *n*. However, higher-order terms need to be included in equations (4), (5) as the number of spheres (or thickness of the sample) increases. This result is similar to other multiple scattering approaches. For example, in the case of a two-beam configuration, it was shown (Cowley & Moodie, 1957) that the scattered beam could be expressed using expansion (17), considering only the scattering terms of the primary beam . The famous two-beam analytical scattering expression is obtained as the sum of the infinite series.

Fig. 9(*a*) shows the scattering probabilities of *n*-times scattering for *ka* = 7, *n*_{p} = 1.01, *ka* = 3*kd* with increasing *N*. The dynamic scattering probability and the kinematic as well as the total scattering probabilities are shown. It seems that all multiple scattering terms increase with the number of spheres, in contrast to the ballistic picture of the (Egerton, 2011), also shown as dashed–dotted lines calculated using the average scattering The main reason for these trends is that the ballistic picture ignores the effect of interference, which drastically affects the overall scattering process.

Fig. 9(*b*) shows the scattering amplitudes *f*_{n}, associated with *n*-times scattering for *ka* = 7, *n*_{p} = 1.01 and *N* = 15. The scattering amplitudes have similar shapes, which mostly contribute to constructive interference, although slight peak misalignments can be seen. The reason for the large-angle scattering, visible on the figure, is due to the small radius used in this example. As mentioned above, a much larger normalized radius would be required to properly account for the low-angle representation of the far-field diffraction pattern. However, because of the comparable scattering of the chosen model and a real carbon atom, the extent of multiple scattering should have been reasonably estimated.

Note that the examples presented here have been performed for electron incident energies of 200 keV, which is a typical value in (*a*) to the left. As a result, spherical harmonics of lower orders would be sufficient to achieve good accuracy. On the other hand, the approximation would lose accuracy at smaller thicknesses, and then full inversion of the system (6) might be required. For higher energies, the forward approximation becomes more accurate but the size of the system increases, making the problem computationally harder.

### 5. Conclusion and perspective

An alternative approach, based on the T-matrix formalism, has been applied to the scattering of fast electrons, by light-atom structures. The validity of important approximations, used in the MS approach, has been discussed, and a multiple scattering approximation framework has been proposed and put in perspective compared with other existing interpretations. While it is worth mentioning that the T-matrix does not scale favourably with large increasing values of the wavelength-normalized radii *ka*, as commonly used for high-energy electron diffraction, it can still provide valuable insights in understanding effects thanks to its ability to solve the wave equation exactly. For a 200 keV electron beam, each spherical coefficient *a*_{p}, *b*_{p} accounts for one wavelength of physical space comparable with the sampling resolution used in typical MS simulations. For a truncation longitudinal number *n* = 400, the total number of coefficients for all atoms would be 2*n*^{2}*N* since azimuthal coefficients are necessary for nonlinear array arrangements. For *N* = 10 000 atoms, this is on the order of a few gigabytes of memory. While inverting a matrix of this size would be computationally intractable, the forward scheme approximation, presented in this paper, makes such a computation viable. Indeed, it offers the possibility of parallel computation for all atoms in the same slice, similar to the way it is traditionally done in finite difference schemes.

In the presented approach, the spherically symmetric effective potential does not model the real electrostatic potential very accurately. However, we anticipate that the proposed interpretation of multiple scattering is equally applicable to the more accurate case of a screened Coulomb potential. Such a potential could naturally be included by using a basis of radial functions representing solutions to the homogeneous Schrödinger's equation, which can be obtained numerically. In fact, since the T-matrix approach mainly relies on the spherical symmetry of the individual scatterers, it can be adapted to any family of radial functions, provided that translational coefficients (15), (16) can be computed numerically.

For periodic structures, the Green's functions (GFs) could be used as basis functions. All GFs are solutions to Schrodinger's equation in a periodic potential, expressed as an expansion upon the lattice vectors of the crystal under consideration. They would automatically satisfy the boundary conditions imposed by the Bloch theorem in crystals. This is the starting point of the Kohn–Rostoker method (Kohn & Rostoker, 1954).

An apparent limitation of both T-matrix and MS approaches lies in the use of the IAM, which, by definition, ignores the effect of bonding. Theoretically, such bonding could be included (Natoli *et al.*, 1986). In fact, the MS approach can be adapted, at a computational cost, to bonding models such as the transferable aspherical atom model (TAAM) or even potentials determined by density functional theory. On the other hand, the single-site T-matrix approach, presented in this paper, cannot work directly with such potentials because the spheres are treated as non-overlapping. A possible option would be to employ effective potentials calculated in the framework of self-consistent electronic structures, where the effect of bonding manifests itself in the bonding-corrected effective potential. However, it is still an open question whether the effects of such bonding play an important role in of organic structures by electron diffraction.

The advantage of the multiple scattering approximation is that it may offer the possibility to include consistently the effect of incoherent and

which can provide a greater insight into effects in real experiments. could, for example, be included in a stochastic fashion as randomly affecting the phase of the scattered waves. might require the use of dedicated electrostatic absorption potentials. It is anticipated that incoherent and may have a significant effect in even when energy filters are used.### 6. Data availability statement

Computer codes, developed in this study, are freely available under the GNU GPL licence v3.0 from https://github.com/ronandrevon/pyScatSpheres. No new ED data were produced in this study.

### APPENDIX A

### T-matrix

The unknown coefficients , are found by imposing the continuity of the wavefunction and its radial derivative at the surface of each sphere *p*:

where is the incident electron wavefunction, the scattered wavefunction of sphere *p* expressed inside the sphere, the scattered wavefunction of sphere *p* expressed outside sphere *p*, and *a*_{p} the radius of sphere *p*.

Using the orthogonality of the spherical harmonics, the following linear system yields the unknown coefficients:

where the translational addition theorem (Dufva *et al.*, 2008) has been used to express the field, scattered by sphere *q* in the reference frame of sphere *p*, formally written as . This operation involves the coupling coefficients , where .

The coefficients *c*_{lm} are related to the incident wave. In the case of a plane wave , , the addition theorem is used to expand the plane wave on the spherical Bessel solutions basis set:

where are the spherical coordinates of the centre of sphere *p* in the global coordinate system and is the angle of incidence with respect to the axis. For propagation in the (*y*,*z*) plane, assumed in this paper, and *e*_{p} is the phase offset at sphere *p*. The different notations are illustrated in Fig. 1.

The coefficients and are expressed as

where .

Each of equations (15), (16) represent *n*^{2} equations for individual spheres where *n* = *l*_{max} is the truncation Note that each equation involves the outside coefficients of all the other spheres through the translational coefficients. As a result, equations (15), (16) represent a *N*×2×*n*^{2} system with as many unknowns as the number of equations and, therefore, can be written in the matrix form (6).

### APPENDIX B

### Multiple scattering picture of the MS approach

The expression for the scattering amplitude *f*(*h*,*k*) of beam *h*,*k* is proportional to the following expression:

where , , , , , is the total thickness of the sample, *N* the number of slices of thickness , *Q*_{h,k,l} = + , *F*_{h,k,l} is the the interaction parameter, ζ the excitation error of beam (*h*,*k*,*l*) and is the excitation error of beam . The excitation error is expressed as

where *a*, *b* and *c* are the lattice constants of the crystal and the wavenumber. Equation (18) expresses the longitudinal distance of beam (*h*,*k*,*l*) to the paraboloid shown in Fig. 10(*a*). This paraboloid is a very accurate representation of the for large wavenumber *K*. Therefore, ζ is very close to the excitation error commonly defined in crystallography.

From equation (17) it is possible to gather terms in powers of into *f*_{h,k}^{(n)}, which corresponds to the incident beam scattered *n* times before contributing to reflections *h*,*k*. The term *Q*_{0,0,0}^{N} only appears for *h* = *k* = 0 and corresponds to the unscattered direct beam. There are *N* terms involving the factor depending on which slice the single scattering event took place. There are terms involving the factor depending on which pair of slices the two-level scattering process is considered, and so on for the higher-order multiple scattering terms (Cowley & Moodie, 1957). This multiple scattering picture is illustrated in 2D in Fig. 10(*b*) for a case with *N* = 3 and using only beams in the zero-order Laue zone (ZOLZ) *l*_{i} = 0.

Using only ZOLZ beams, *f*_{h,k}^{(1)}, is expanded as

which is the well known kinematic scattering regime, where the sum converges to the standard

curvature factor for infinitely thick slices , , .For a two-level scattering using only ZOLZ beams, *f*_{h}^{(2)} would expand as (written in 2D here)

### Acknowledgements

This research was conducted at the UK Research and Innovation (UKRI) Research Complex at Harwell, Science and Technology Facilities Council (STFC) Rutherford Appleton Laboratory.

### Funding information

This work was supported by the Biotechnology and Biological Sciences Research Council (BBSRC) grant No. BB/S007040/1.

### References

Allen, L. J., D'Alfonso, A. J. & Findlay, S. D. (2015). *Ultramicroscopy*, **151**, 11–22. Web of Science CrossRef CAS PubMed Google Scholar

Balanis, C. (1989). *Advanced Engineering Electromagnetics.* New York: John Wiley and Sons. Google Scholar

Bethe, H. A. (1928). *Ann. Phys.* **392**, 55–129. CrossRef Google Scholar

Chen, J. H., Van Dyck, D. & Op de Beeck, M. (1997). *Acta Cryst.* A**53**, 576–589. CrossRef CAS Web of Science IUCr Journals Google Scholar

Chuang, S. L. (1996). *Physics of Optoelectronic Devices*, edited by J. W. Goodman, pp. 631–650, Wiley Series in Pure and Applied Optics. New York: John Wiley and Sons. Google Scholar

Cowley, J. M. (1995). *Diffraction Physics*, 3rd ed. Amsterdam: North Holland Publishing Company. Google Scholar

Cowley, J. M. & Moodie, A. F. (1957). *Acta Cryst.* **10**, 609–619. CrossRef IUCr Journals Web of Science Google Scholar

Dederichs, P. H. (1971). *Berichte der Kernforschungsanlage* Jülich No. 797 (printed as manuscript). Google Scholar

Drevon, T. R. (2021). *pyScatSpheres*. https://pypi.org/project/pyScatSpheres/. Google Scholar

Dufva, T. J., Sarvas, J. & Sten, J. C. (2008). *Prog. Electromagn. Res. B*, ** 4**, 79–99. CrossRef Google Scholar

Egel, A., Pattelli, L., Mazzamuto, G., Wiersma, D. S. & Lemmer, U. (2017). *J. Quant. Spectrosc. Radiat. Transfer*, **199**, 103–110. Web of Science CrossRef CAS Google Scholar

Egerton, R. (2011). In *Electron Energy-loss Spectroscopy in the Electron Microscope*. New York: Springer. Google Scholar

Eremin, J. A., Orlov, N. V. & Rozenberg, V. I. (1995). *J. Atmos. Terrestrial Phys.* **57**, 311–319. CrossRef Web of Science Google Scholar

Gemmi, M., Mugnaioli, E., Gorelik, T. E., Kolb, U., Palatinus, L., Boullay, P., Hovmoller, S. & Abrahams, J. P. (2019). *ACS Cent. Sci.* **5**, 1315–1329. Web of Science CrossRef CAS PubMed Google Scholar

Glaeser, R. M. & Downing, K. H. (1993). *Ultramicroscopy*, **52**, 478–486. CrossRef CAS PubMed Web of Science Google Scholar

Godin, O. A. (2011). *J. Acoust. Soc. Am.* **130**, EL135–EL141. Web of Science CrossRef PubMed Google Scholar

Hamid, A.-K., Ciric, I. R. & Hamid, M. (1990*a*). *Can. J. Phys.* **68**, 1419–1428. CrossRef Web of Science Google Scholar

Hamid, A.-K., Ciric, I. R. & Hamid, M. (1990*b*). *Can. J. Phys.* **68**, 1157–1165. CrossRef Web of Science Google Scholar

Hirsch, P. B., Howie, A., Nicholson, R. B., Pashley, D. W. & Whelan, M. J. (1965). *Electron Microscopy of Thin Crystals*. London: Butterworths. Google Scholar

Howie, A. (1963). *Proc. R. Soc. Lond. A*, **271**, 268–287. CAS Google Scholar

Howie, A. (2014). *J. Phys. Conf. Ser.* **522**, 012001. Google Scholar

Ishizuka, K. (1982). *Acta Cryst.* A**38**, 773–779. CrossRef CAS Web of Science IUCr Journals Google Scholar

Ishizuka, K. & Uyeda, N. (1977). *Acta Cryst.* A**33**, 740–749. CrossRef CAS IUCr Journals Web of Science Google Scholar

Kirkland, E. J. (2019). *Advanced Computing in Electron Microscopy*, 3rd ed. Cham: Springer. Google Scholar

Klar, P. B., Xu, H., Steciuk, G., Cho, J., Zou, X., Palatinus, L. & Republic, C. (2021). *chemrxiv*, pp. 1–25. Google Scholar

Kohn, W. & Rostoker, N. (1954). *Phys. Rev.* **94**, 1111–1120. CrossRef CAS Web of Science Google Scholar

Korringa, J. (1947). *Physica*, **13**, 392–400. CrossRef Web of Science Google Scholar

Korringa, J. (1994). *Phys. Rep.* **238**, 341–360. CrossRef Web of Science Google Scholar

Kottke, A. (2020). *pygmm*. https://zenodo.org/badge/latestdoi/21452/arkottke/pygmm. Google Scholar

Latychevskaia, T. & Abrahams, J. P. (2019). *Acta Cryst.* B**75**, 523–531. Web of Science CrossRef IUCr Journals Google Scholar

Metherell, A. & Fisher, R. (1969). *Phys. Status Solidi B*, **32**, 551–562. CrossRef CAS Web of Science Google Scholar

Moine, O. & Stout, B. (2005). *J. Opt. Soc. Am. B*, **22**, 1620. Web of Science CrossRef Google Scholar

Müller, H. J. (1965). *Physica*, **31**, 688–692. Google Scholar

Nannenga, B. L. & Gonen, T. (2019). *Nat. Methods*, **16**, 369–379. Web of Science CrossRef CAS PubMed Google Scholar

Nannenga, B. L., Shi, D., Leslie, A. G. W. & Gonen, T. (2014). *Nat. Methods*, **11**, 927–930. Web of Science CrossRef CAS PubMed Google Scholar

Natoli, C. R., Benfatto, M. & Doniach, S. (1986). *Phys. Rev. A*, **34**, 4682–4694. CrossRef CAS PubMed Web of Science Google Scholar

Oleynikov, P., Hovmöller, S. & Zou, X. D. (2007). *Ultramicroscopy*, **107**, 523–533. Web of Science CrossRef PubMed CAS Google Scholar

Ophus, C. (2017). *Adv. Struct. Chem. Imaging*, **3**, 1–11. CrossRef PubMed Google Scholar

Palatinus, L., Jacob, D., Cuvillier, P., Klementová, M., Sinkler, W. & Marks, L. D. (2013). *Acta Cryst.* A**69**, 171–188. Web of Science CrossRef CAS IUCr Journals Google Scholar

Rehr, J. J., Kas, J. J., Prange, M. P., Sorini, A. P., Takimoto, Y. & Vila, F. (2009). *C. R. Phys.* **10**, 548–559. Web of Science CrossRef CAS Google Scholar

Rehr, J. J., Kas, J. J., Vila, F. D., Prange, M. P. & Jorissen, K. (2010). *Phys. Chem. Chem. Phys.* **12**, 5503–5513. Web of Science CrossRef CAS PubMed Google Scholar

Schiff, L. I. (1955). *Quantum Mechanics*, 2nd ed. New York: McGraw-Hill. Google Scholar

Sébilleau, D., Gunnella, R., Wu, Z. Y., Matteo, S. D. & Natoli, C. R. (2006). *J. Phys. Condens. Matter*, **18**, R175–R230. Google Scholar

Sébilleau, D., Hatada, K. & Ebert, H. (2017). *Multiple Scattering Theory for Spectroscopies.* Cham: Springer. Google Scholar

Sébilleau, D., Natoli, C., Gavaza, G. M., Zhao, H., Da Pieve, F. & Hatada, K. (2011). *Comput. Phys. Commun.* **182**, 2567–2579. Google Scholar

Silva, G. T., Baggio, A. L., Lopes, J. H. & Mitri, F. G. (2012). arXiv:1210.2116. Google Scholar

Subramanian, G., Basu, S., Liu, H., Zuo, J. & Spence, J. C. H. (2015). *Ultramicroscopy*, **148**, 87–93. Web of Science CrossRef CAS PubMed Google Scholar

Vainshtein, B. (1964). *Structure Analysis by Electron Diffraction*, pp. 114–204. Oxford: Pergamon Press. Google Scholar

Van Dyck, D. & Op de Beeck, M. (1996). *Ultramicroscopy*, **64**, 99–107. CrossRef CAS Web of Science Google Scholar

Zabloudil, J., Hammerling, R., Weinberger, P. & Szunyogh, L. (2005). Editors. *Electron Scattering in Solid Matter*. Berlin: Springer. Google Scholar

Zuo, J. M. & Weickenmeier, A. L. (1995). *Ultramicroscopy*, **57**, 375–383. CrossRef CAS Web of Science Google Scholar

This is an open-access article distributed under the terms of the Creative Commons Attribution (CC-BY) Licence, which permits unrestricted use, distribution, and reproduction in any medium, provided the original authors and source are cited.