research papers\(\def\hfill{\hskip 5em}\def\hfil{\hskip 3em}\def\eqno#1{\hfil {#1}}\)

Journal logoFOUNDATIONS
ADVANCES
ISSN: 2053-2733

Dynamical diffraction of high-energy electrons by light-atom structures: a multiple forward scattering interpretation

crossmark logo

aSTFC, Rutherford Appleton Laboratory, Didcot, OX11 0FA, United Kingdom, and bCCP4, Research Complex at Harwell, Rutherford Appleton Laboratory, Didcot, OX11 0FA, United Kingdom
*Correspondence e-mail: eugene.krissinel@stfc.ac.uk

Edited by L. Palatinus, Czech Academy of Sciences, Czech Republic (Received 21 March 2022; accepted 8 December 2022; online 9 February 2023)

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 dynamical diffraction 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 forward scattering 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.

1. Introduction

1.1. Motivation

The theory of dynamical diffraction was first developed in 1928 (Bethe, 1928[Bethe, H. A. (1928). Ann. Phys. 392, 55-129.]; Cowley, 1995[Cowley, J. M. (1995). Diffraction Physics, 3rd ed. Amsterdam: North Holland Publishing Company.]), 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[Korringa, J. (1947). Physica, 13, 392-400.], 1994[Korringa, J. (1994). Phys. Rep. 238, 341-360.]; Kohn & Rostoker, 1954[Kohn, W. & Rostoker, N. (1954). Phys. Rev. 94, 1111-1120.]; Dederichs, 1971[Dederichs, P. H. (1971). Berichte der Kernforschungsanlage Jülich No. 797 (printed as manuscript).]) has proven to be very efficient for describing electronic properties of matter. Besides, MST provides an intuitive picture of dynamical diffraction. Partial wave scattering theory (Schiff, 1955[Schiff, L. I. (1955). Quantum Mechanics, 2nd ed. New York: McGraw-Hill.]) provides another description of multiple scattering with the possibility to take into account inelastic scattering (Howie, 1963[Howie, A. (1963). Proc. R. Soc. Lond. A, 271, 268-287.]) and loss of coherency (Howie, 2014[Howie, A. (2014). J. Phys. Conf. Ser. 522, 012001. ]). Other less commonly known attempts at an intuitive depiction of dynamical diffraction have been proposed to explain electron channelling in crystals, for example by VanDyck & Op de Beeck (1996[Van Dyck, D. & Op de Beeck, M. (1996). Ultramicroscopy, 64, 99-107. ]). Dynamical diffraction 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[Hirsch, P. B., Howie, A., Nicholson, R. B., Pashley, D. W. & Whelan, M. J. (1965). Electron Microscopy of Thin Crystals. London: Butterworths.]; Glaeser & Downing, 1993[Glaeser, R. M. & Downing, K. H. (1993). Ultramicroscopy, 52, 478-486.]; Subramanian et al., 2015[Subramanian, G., Basu, S., Liu, H., Zuo, J. & Spence, J. C. H. (2015). Ultramicroscopy, 148, 87-93.]). 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 structure determination process. However, structures have been successfully determined from ED (electron diffraction) data using the standard kinematic theory (Gemmi et al., 2019[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.]; Nannenga et al., 2014[Nannenga, B. L., Shi, D., Leslie, A. G. W. & Gonen, T. (2014). Nat. Methods, 11, 927-930.]; Nannenga & Gonen, 2019[Nannenga, B. L. & Gonen, T. (2019). Nat. Methods, 16, 369-379. ]). Although dynamical refinement (Palatinus et al., 2013[Palatinus, L., Jacob, D., Cuvillier, P., Klementová, M., Sinkler, W. & Marks, L. D. (2013). Acta Cryst. A69, 171-188.]) usually leads to better intensity predictions (Gemmi et al., 2019[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.]), the agreement between theory and experiment is still significantly worse than that obtained for X-ray data (Oleynikov et al., 2007[Oleynikov, P., Hovmöller, S. & Zou, X. D. (2007). Ultramicroscopy, 107, 523-533. ]; Klar et al., 2021[Klar, P. B., Xu, H., Steciuk, G., Cho, J., Zou, X., Palatinus, L. & Republic, C. (2021). chemrxiv, pp. 1-25.]).

Although there is no quantitative prediction of the effects of dynamical diffraction on the ability to solve a crystal structure, it may be speculated that dynamical diffraction 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[Sébilleau, D., Hatada, K. & Ebert, H. (2017). Multiple Scattering Theory for Spectroscopies. Cham: Springer.]; Zabloudil et al., 2005[Zabloudil, J., Hammerling, R., Weinberger, P. & Szunyogh, L. (2005). Editors. Electron Scattering in Solid Matter. Berlin: Springer.]). However, in the context of high-energy electron diffraction, the multislice (MS) (Cowley & Moodie, 1957[Cowley, J. M. & Moodie, A. F. (1957). Acta Cryst. 10, 609-619.]) approach and Bloch-wave (BW) (Bethe, 1928[Bethe, H. A. (1928). Ann. Phys. 392, 55-129.]; Metherell & Fisher, 1969[Metherell, A. & Fisher, R. (1969). Phys. Status Solidi B, 32, 551-562.]) 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[Ishizuka, K. & Uyeda, N. (1977). Acta Cryst. A33, 740-749.]). 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[Ishizuka, K. (1982). Acta Cryst. A38, 773-779.]; Chen et al., 1997[Chen, J. H., Van Dyck, D. & Op de Beeck, M. (1997). Acta Cryst. A53, 576-589.]). 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[Kirkland, E. J. (2019). Advanced Computing in Electron Microscopy, 3rd ed. Cham: Springer.]). A few very efficient MS implementations have been reported to date (Ophus, 2017[Ophus, C. (2017). Adv. Struct. Chem. Imaging, 3, 1-11.]). Some of them are targeted at convergent-beam electron diffraction (CBED), while others have special features such as inclusion of inelastic scattering (Allen et al., 2015[Allen, L. J., D'Alfonso, A. J. & Findlay, S. D. (2015). Ultramicroscopy, 151, 11-22.]). 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 structure determination 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[Zuo, J. M. & Weickenmeier, A. L. (1995). Ultramicroscopy, 57, 375-383.]) have been proposed to include a large number of beams, BW can hardly be applied to very large structures since the scattering matrix becomes prohibitively large. Due to the unfavourable O(N3) 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., 1990a[Hamid, A.-K., Ciric, I. R. & Hamid, M. (1990a). Can. J. Phys. 68, 1419-1428.],b[Hamid, A.-K., Ciric, I. R. & Hamid, M. (1990b). Can. J. Phys. 68, 1157-1165.]; Eremin et al., 1995[Eremin, J. A., Orlov, N. V. & Rozenberg, V. I. (1995). J. Atmos. Terrestrial Phys. 57, 311-319.]), optics (Moine & Stout, 2005[Moine, O. & Stout, B. (2005). J. Opt. Soc. Am. B, 22, 1620.]) and acoustics (Silva et al., 2012[Silva, G. T., Baggio, A. L., Lopes, J. H. & Mitri, F. G. (2012). arXiv:1210.2116.]; Godin, 2011[Godin, O. A. (2011). J. Acoust. Soc. Am. 130, EL135-EL141.]). 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, abinitio 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[Rehr, J. J., Kas, J. J., Prange, M. P., Sorini, A. P., Takimoto, Y. & Vila, F. (2009). C. R. Phys. 10, 548-559. ]). Similar approaches based on the T-matrix also exist in electron energy-loss spectroscopy (Sébilleau et al., 2006[Sébilleau, D., Gunnella, R., Wu, Z. Y., Matteo, S. D. & Natoli, C. R. (2006). J. Phys. Condens. Matter, 18, R175-R230. ]).

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 forward scattering approximation. We propose a comparison with the MS interpretation of multiple scattering. In addition, we discuss the validity of the forward scattering 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[Kirkland, E. J. (2019). Advanced Computing in Electron Microscopy, 3rd ed. Cham: Springer.]):

[\left[-{{\hbar^{2}} \over {2m_{\rm e}}}\boldnabla_{\bf r}^{2}-eV({\bf r})\right]\Psi = E\Psi, \eqno(1)]

where [\hbar] is Planck's constant, [m_{\rm e}] the mass of the electron, e the elementary charge, [V({\bf r})] the spatially varying electrostatic potential created by the atoms, and E the electron energy.

In an ED experiment, [E\gg V(r\to\infty) = 0], which corresponds to the continuum, non-quantized state of the system (Chuang, 1996[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.]). Therefore, we solve (1)[link] with E as a parameter, equal to the energy of incident electrons.

The total wavefunction Ψ is described as the sum of an incident wave [\Psi^{({\rm i})}] and a scattered wave [\Psi^{({\rm s})}]:

[\Psi({\bf r}) = \Psi^{({\rm i})}+\Psi^{({\rm s})}.\eqno(2)]

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 [\Psi^{({\rm s})}] 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 [k_{0} = 2\pi/\lambda] (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 V0 should provide similar values of scattering cross sections. The scattering cross section 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[link]. This formulation is well established (Hamid et al., 1990a[Hamid, A.-K., Ciric, I. R. & Hamid, M. (1990a). Can. J. Phys. 68, 1419-1428.],b[Hamid, A.-K., Ciric, I. R. & Hamid, M. (1990b). Can. J. Phys. 68, 1157-1165.]; Eremin et al., 1995[Eremin, J. A., Orlov, N. V. & Rozenberg, V. I. (1995). J. Atmos. Terrestrial Phys. 57, 311-319.]) 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:

[\eqalign{\boldnabla_{{\bf r}_{p}}^{2}\Psi+k_{p}^{2}\Psi &= 0\cr k_{p} &= k_{0}n_{p}\cr k_{0}& = \sqrt{{{2m_{\rm e}E} \over {\hbar^{2}}} }\cr n_{p}& = \sqrt{1+{{V_{p}} \over {E}} },}]

where [V_{p}\geq 0] is the constant positive potential inside sphere p of radius ap centred at [{\bf d}_{p}], kp the wavenumber inside the sphere and np is sometimes referred to as the refractive index by analogy with optics. In the remainder, np will be referred to as the potential strength for clarity purposes. The scattered wavefunction [\Psi^{({\rm s})}] can be decomposed into the scattered wavefunctions of individual atom spheres,

[\Psi^{({\rm s})} = \textstyle\sum\limits_{p = 1}^{N}\Psi_{p}^{({\rm s})},\eqno(3)]

where N is the number of spheres. The p-sphere scattered wave [\Psi_{p}^{({\rm s})}] can itself be partitioned into a part inside sphere p noted [\Psi_{p}^{({\rm in})}] and a part outside sphere p noted [\Psi_{p}^{({\rm out})}]. The expressions for [\Psi_{p}^{({\rm in})}] and [\Psi_{p}^{({\rm out})}] are detailed as

[\Psi_{p}^{({\rm in})}({\bf r}_{p}) = \textstyle\sum\limits_{l = 0}^{\infty}j_{l}(k_{p}r_{p}) \sum\limits_{m = -l}^{m = l}a_{p\semi lm}Y_{l}^{m}(\theta_{p},\phi_{p})\eqno(4)]

[\Psi_{p}^{({\rm out})}({\bf r}_{p}) = \textstyle\sum\limits_{l = 0}^{\infty}h_l^{(1)}(k_{0}r_{p}) \sum\limits_{m = -l}^{m = l}b_{p\semi lm}Y_{l}^{m}(\theta_{p},\phi_{p}),\eqno(5)]

where [p = 1\ldots N], jl and hl(1) are the spherical Bessel and Hankel functions of the first kind, respectively, Ylm 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 [{\bf r}_{p}].

[Figure 1]
Figure 1
Setup and conventions for the scattering problem solved with the T-matrix approach. A plane wave with a wavevector k0 with incident angles ([\theta_{\rm i}], [\phi_{\rm i}]) interacts with N spherical atoms, each with individual location (for example, atom p at dp), radius ( ap) and electrostatic potential strength ( np). The scattered wave is observed at observation point P at r in the global reference frame and at individual distances in the local reference frame of the spheres (for example, [|{\bf r}_{p}|]). See discussion in the text for additional details. In practice, the incident plane wave is taken along the +z axis and the structure oriented accordingly.

The unknown coefficients [a_{p\semi lm}], [b_{p\semi lm}] 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[link]), this results in the following linear system of equations:

[({\bf I}-{\bf T}){\bf A} = {\bf L},\eqno(6)]

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)[link].

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 [|f(\theta,\phi)|^{2}] 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 [f(\theta,\phi)] no longer depends on the distance to the scattering specimen. Using the asymptotic behaviour hl(1)(k0rp) [\simeq (-j)^{l+1}\{[\exp(jk_{0}r_{p})] / (k_{0}r_{p})\}] and since [\theta_{p}\simeq\theta], [\phi_{p}\simeq\phi], the far-field scattering amplitude [f_{p}(\theta,\phi)] from each individual sphere p can be written as

[f_{p}(\theta,\phi) = \textstyle\sum\limits_{l = 0}^{\infty}\sum\limits_{m = -l}^{l}(-j)^{l+1}b_{p\semi lm}Y_{l}^{ m}(\theta,\phi).\eqno(7)]

The total far-field scattering amplitude [f(\theta,\phi)] 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 ap and bp. Since, in the far field, [r_{p}\simeq r-{\bf d}_{p}\cdot{\bf e}_{r}],

[f(\theta,\phi) = \textstyle\sum\limits_{p = 1}^{N}f_{p}(\theta,\phi)\exp(-jk_{0}{\bf e}_{r}\cdot{\bf d}_{ p}).\eqno(8)]

The normalized differential scattering cross section can be obtained as

[{{\sigma(\theta,\phi)} \over {\pi a_{p}^{2}}} = {{4\pi r^{2}} \over {\pi a_{p}^{2}}} \Big\|{{{\Psi^{({\rm s})}(r,\theta,\phi)} \over {\Psi^{({\rm i})}(r,\theta,\phi)}}}\Big\|^{2} = {{ 4|f(\theta,\phi)|^{2}} \over {\left(k_{0}a_{p}\right)^{2}}},\eqno(9)]

where we have used [\Psi^{({\rm s})}(r,\theta,\phi)] [\simeq_{r\rightarrow\infty}[\exp(jk_{0} r)/ (k_{0}r)]f(\theta,\phi)].

3. Real-space forward multiple scattering picture

3.1. T-matrix forward scattering approximations

Equation (6)[link] 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 [a_{p\semi lm} = c_{lm}u_{p\semi l}], [b_{p\semi lm} = c_{lm}v_{p\semi l}]. These are indeed the well known analytical solutions of Mie scattering by a soft sphere (Balanis, 1989[Balanis, C. (1989). Advanced Engineering Electromagnetics. New York: John Wiley and Sons.]).

The cross-coupling matrix T accounts for multiple scattering effects. If A is written as [{\bf A} = (\ldots a_{p\semi lm},b_{p\semi lm}\ldots)^{\rm T}], then T is a matrix with only off-diagonal components:

[{\bf T} = \left[\matrix{{\bf 0}&\ldots{\bf T}_{1q}\ldots &\ldots{\bf T}_{1p}\ldots &{\bf T}_{1N}\cr {\bf T}_{q1}&\ldots{\bf 0}\ldots &\ldots{\bf T}_{qp}\ldots &{\bf T}_{qN}\cr {\bf T}_{p1}&\ldots{\bf T}_{pq}\ldots &\ldots{\bf 0}\ldots &{\bf T}_{pN}\cr {\bf T}_{N1}&\ldots{\bf T}_{Nq}\ldots &\ldots {\bf T}_{Np}\ldots &{\bf 0}\cr }\right],\eqno(10)]

where Tpq represents the scattering from sphere p due to the scattering from sphere q. By considering all components [{\bf T}_{pq}], [{\bf T}_{qp}], we assume that the scattering from sphere p affects scattering from sphere q and vice versa. Matrix (10)[link] 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)[link] requires matrix inversion which is computationally expensive.

However, as detailed below in Section 4.2[link], backscattering can be neglected for very fast electrons, which is known as the forward scattering approximation. This results in T being lower triangular if the spheres are sorted in ascending order along [{\bf e}_{z}]. 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 forward scattering 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 forward scattering approximation converts the implicit self-consistent T-matrix scheme into a forward scattering explicit scheme similar to the approach taken in the MS method.

3.2. Multiple scattering approximations

Since [{\bf A}_{0} = {\bf L}] represents single scattering, we can establish that [{\bf A}_{1} = {\bf TA}_{0}] accounts for secondary scattering. Similarly, outward scattering amplitudes from sphere p can be written as

[\eqalignno{{\bf b}_{p} &= {\bf b}_{p}^{(0)}+\textstyle\sum\limits_{q,z_{q}\lt z_{p}}{\bf T}_{pq}{\bf b}_{q} ^{(0)}&\cr &\quad +{\bf T}_{pq}\textstyle\sum\limits_{q,z_{q}\lt z_{p}}{\bf T}_{qr}\sum\limits_{r,z_{r}\lt z_{q}}{\bf b}_{r }^{(0)}+\ldots &\cr &= \textstyle\sum\limits_{n = 0}{\bf T}^{n}{\bf b}_{p}^{(0)} = \sum\limits_{n = 1}{\bf b}_{p}^{(n)},&(11)}]

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[Korringa, J. (1947). Physica, 13, 392-400.], 1994[Korringa, J. (1994). Phys. Rep. 238, 341-360.]; Kohn & Rostoker, 1954[Kohn, W. & Rostoker, N. (1954). Phys. Rev. 94, 1111-1120.]). This forward multiple scattering picture is illustrated in Fig. 2[link] for the case of a two-atom system. From a computational point of view, equation (11)[link] may offer an advantage over the full forward scattering approximation if only a few n-time scattering terms are necessary. Indeed, since each term in the expansion depends on precomputed coefficients [b_{p\semi lm}^{(0)}], computation of the [b_{p\semi lm}^{(n)}] coefficients for the nth scattering term can be massively parallelized for each atom.

[Figure 2]
Figure 2
Two-level forward scattering approximation using the T-matrix approach with an incident plane wave propagating upwards. The third column represents scattering from the top scatterer, located at distance kd from the bottom scatterer. By neglecting backscattering, this can be approximated by the sum of (i) the kinematic scattering of the incident plane wave by the top scatterer (first column) and (ii) the secondary scattering following kinematic scattering from the bottom scatterer (second column). The scattering from the bottom scatterer is viewed by the top scatterer as a particular incident wave (neither planar nor spherical) emitted from a point source at a distance kd.

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[Egerton, R. (2011). In Electron Energy-loss Spectroscopy in the Electron Microscope. New York: Springer.]) that the probability of an electron being elastically scattered n times after passing through an amorphous sample of thickness z follows a Poisson distribution [P_{\rm e}^{n} = ({1} / {n!})(z/l_{\rm e})^{n}\exp(-z/l_{\rm e})], where [l_{\rm e}] is the elastic mean free path. It is paramount to note that this classical expression models incoherent scattering, 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 incoherent scattering if inelastic scattering is strongly present. Indeed, inelastic scattering not only induces a transfer of energy between scatterers but also a loss of coherency through phase shifts (Howie, 2014[Howie, A. (2014). J. Phys. Conf. Ser. 522, 012001. ]). As a result, quantum systems, where inelastic scattering is predominant, would scatter incoherently. For example, inelastic scattering 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 elastic scattering 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 cross section as follows. Let 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 [|A|^{2} = 1/Sz] 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[Vainshtein, B. (1964). Structure Analysis by Electron Diffraction, pp. 114-204. Oxford: Pergamon Press.]) is defined as [J_{0} = \hbar k_{0}/m|A|^{2} = v_{0}/Sz], where v0 is the speed of the incident electrons.

By definition of the cross section, [J_{0}\sigma] is the probability of an electron being scattered every second. Since it takes a time [\Delta t = z/v_{0}] for an electron to pass through the specimen, the overall probability of a single electron being scattered is therefore [P_{\rm scat} = \sigma J_{0}\Delta t]. Therefore, after inserting the expressions for J0 and [\Delta t],

[P_{\rm scat} = {{\sigma} \over {S}},]

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

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

[\eqalignno{ f(\theta)& = \textstyle\sum\limits_{n = 1}^{N}f_{n}(\theta)&\cr \sigma &= \textstyle\sum\limits_{n = 1}^{N}\sigma_{n}&\cr \sigma_{n}& = | f_{n}|^{2}+2{\rm Re} \left(f_{n}\textstyle\sum\limits_{m\gt n}^{N} f_{m}^{*}\right).&(12)}]

We can then define [P_{n} = \sigma_{n}/S] as the probability of an electron being scattered n times. The choice of n-time scattering cross section [\sigma_{n}] is arbitrary precisely because of the cross-coupling interference terms fnfm*. 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[Cowley, J. M. & Moodie, A. F. (1957). Acta Cryst. 10, 609-619.]). Our definition of scattering probabilities should reflect this behaviour.

The probability of an electron not being scattered is naturally [P_{\rm coh} = 1-P_{\rm scat}], where [P_{\rm scat} = \sum_{n}P_{n}]. Note that, since [S\gg\sigma], 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 cross section [\sigma_{\rm a} = \sigma/N]. The corresponding scattering probability is then P(dz) = [\sigma_{\rm a}\rho dz = dz/l_{\rm e}] where [\rho = 1/Sdz] is the density (since Sdz contains only one scatterer) and [l_{\rm e} = 1/\sigma_{\rm a}\rho] is the mean free path. This is consistent with the definition used by Latychevskaia & Abrahams (2019[Latychevskaia, T. & Abrahams, J. P. (2019). Acta Cryst. B75, 523-531.]).

3.4. First Born approximation and kinematic scattering

Let us assume normal incidence with [\theta_{\rm i} = 0], so that ep = exp(jk0dp), and note that the Ewald sphere can be represented by q = [K(\sin\theta\cos\phi,\sin\theta\sin\phi,1-\cos\theta)] in reciprocal space (or transfer vector momentum space), with [K = 1/\lambda] being the wavenumber in crystallographic convention. Keeping only the kinematic term in equation (11)[link] and inserting in the far field in equation (8)[link] give

[f(\theta,\phi) = \textstyle\sum\limits_{p = 1}^{N}f_{p}^{(e)}\exp(2j\pi{\bf q}\cdot{\bf d}_{p}),\eqno(13)]

where fp(e) is the scattering amplitude of sphere p in the first Born approximation, also known as the form factor. Equation (13)[link] corresponds to the standard kinematic formula for the structure factor 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[link] 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 forward scattering 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[Cowley, J. M. & Moodie, A. F. (1957). Acta Cryst. 10, 609-619.]), which, although analogous to the one presented above, differs in that it is stated in reciprocal space. More details on the differences between the two interpretations can be found in Appendix B[link].

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[Egel, A., Pattelli, L., Mazzamuto, G., Wiersma, D. S. & Lemmer, U. (2017). J. Quant. Spectrosc. Radiat. Transfer, 199, 103-110.]; Kottke, 2020[Kottke, A. (2020). pygmm. https://zenodo.org/badge/latestdoi/21452/arkottke/pygmm.]). There are also packages targeted at photon and electron spectroscopy and related to the T-matrix, such as msSpec (Sébilleau et al., 2011[Sébilleau, D., Natoli, C., Gavaza, G. M., Zhao, H., Da Pieve, F. & Hatada, K. (2011). Comput. Phys. Commun. 182, 2567-2579. ]) and FEFF9 (Rehr et al., 2010[Rehr, J. J., Kas, J. J., Vila, F. D., Prange, M. P. & Jorissen, K. (2010). Phys. Chem. Chem. Phys. 12, 5503-5513.]). Here, an open-source package, adapted to high-energy electron diffraction and including the forward scattering approximation, has been developed and made available (Drevon, 2021[Drevon, T. R. (2021). pyScatSpheres. https://pypi.org/project/pyScatSpheres/.]).

4.1. Validity of the implementation

In practice, the size of the matrix in equation (6)[link] has to be truncated to a maximum integer order [\nu_{\max}] 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 kamax = maxp kap where ap is the radius of the spheres and k the wavenumber.

Besides, the translational addition theorem (Dufva et al., 2008[Dufva, T. J., Sarvas, J. & Sten, J. C. (2008). Prog. Electromagn. Res. B, 4, 79-99.]) 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 [\nu_{\max}]. Once again, the radius of the spheres determines the choice of [\nu_{\max}]. This is illustrated in Fig. 3[link](a), where the error between hl(1) Y42, computed at the origin, using the translational addition theorem with [{\bf d}_{p} = (0,3,5)] and [\nu_{\max} = 10], is displayed in log scale.

[Figure 3]
Figure 3
(a) Error difference in evaluating h2(1)Y42 using (i) the standard product of the radial Hankel function h2(1) and spherical harmonic Y42 and (ii) the translational addition theorem (14)[link], (15)[link] at a distance [{\bf d}_{p} = (0,3,5)] (in normalized radii) from the origin with expansion order [\nu_{\max} = 10]. The circle represents spheres of radius 1. Colour axis in log scale. (b) Continuity error at the sphere boundaries [i.e. difference between the solution computed at the boundary using equations (4)[link] and (5)[link], integrated over the sphere boundary]. For this example, N = 4 and kamax = 4. The error reaches machine accuracy at expansion order [\nu_{\max} = 10].

The choice of [\nu_{\max}] is best evaluated by assessing the continuity of the wavefunction at the surface of the spheres as shown in Fig. 3[link](b) for N = 4 and kamax = 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 forward scattering 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[Kirkland, E. J. (2019). Advanced Computing in Electron Microscopy, 3rd ed. Cham: Springer.]) which results in the potential shown in Fig. 4[link](a).

[Figure 4]
Figure 4
(a) Electrostatic potential, created by an atom of carbon in the IAM (blue solid line). The blue patches show a multi-shell approximation model, which can be used to represent the potential in a T-matrix approach. (b) Scattering amplitude in the Born approximation for the multi-shell model with increasing truncation radius kamax.

The solution to Schrödinger's equation in such a potential can only be solved perturbatively (Müller, 1965[Müller, H. J. (1965). Physica, 31, 688-692.]) and is beyond the scope of this study. However, indicative values for kp and kap can be used with a multi-shell representation as shown in the form of blue patches in Fig. 4[link](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 forward scattering approximation. In Fig. 4[link](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)[link], (5)[link] would need to be included. However, the medium truncation range should be sufficient for providing a reasonable picture of dynamical scattering. Fig. 5[link](a) shows the total scattering cross section of a single sphere with increasing radius for a range of values of potential strength np. The dark blue curve shows the locations of the spherical shell for carbon. It becomes almost flat for the parameter set ( ka = 30, np = 1.001), from which it keeps increasing very slightly to reach the asymptotic value [\sigma\simeq 0.0035]. This value is almost identical to the average elastic cross section for real carbon in the Born approximation (Latychevskaia & Abrahams, 2019[Latychevskaia, T. & Abrahams, J. P. (2019). Acta Cryst. B75, 523-531.]) 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.

[Figure 5]
Figure 5
(a) Total scattering cross section σ at 200 keV for a few values of np over a range of normalized radii ka. The dark blue curve shows the location of the carbon spherical shells. (b) Shape of the scattering amplitude for normalized radius. The green curves correspond to ka = 3 and potential's strength np = 1.1. The blue curves correspond to ka = 10, np = 1.02 which is one of the carbon shells. Dashed lines show the Born approximation and dashed–dotted lines show the phase grating approximation used in MS. The Born approximation is very good for the carbon data point used.

Fig. 5[link](b) shows the shape of the far-field amplitudes for two sets of ka, np values. The first one is for the extreme value of ka = 3 and np = 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[link](a) shows a (ka,np) 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

[{\rm err}(b_{p\semi l,m}) = \textstyle\sum\limits_{p\semi l,m}| b_{p\semi l,m}-b_{p\semi l,m}^{\rm approx}|. \eqno(14)]

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

[Figure 6]
Figure 6
Error of [b_{p\semi lm}] [equation (14)[link]] in (ka,np) space (colour axis in log scale) for a two-scatterer system using the (a) kinematic scattering approximation, (b) forward scattering approximation. The blue dots correspond to the location of the spherical shells of a carbon atom at E = 200 keV (shown only for first five shells for visualization purposes).

For this case, a value of kd = 3ka is used, which is close to interatomic distances. Overall, it is clear that the forward scattering 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 forward scattering 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 np results in an overall good approximation of both the uncoupled and forward scattering 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 cross section. On the other hand, the forward scattering approximation improves with larger values of ka, which make backward scattering less likely.

Fig. 7[link](a) shows the scattering coefficients [b_{p\semi l,m}] in complex space for the given two-body problem, where np = 1.01, kd = 3ka, 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 forward scattering approximation and the kinematic approximation. In order to get this picture, the coefficients were divided by the phase factor at sphere p, i.e. [b_{p\semi lm}\exp(jkd)]. The arms of the spiral are rotated by [\pi/2] from one order to the other as a result of the jl factor in the spherical expansion of the incident plane wave. All the coefficients increase with increasing ka as a result of the increasing scattering cross section 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.

[Figure 7]
Figure 7
(a) Complex space scattering coefficients [b_{p\semi l,m}] for a two-body problem for np = 1.01, kd = 3ka with increasing radius ka. Only the coefficients for sphere p = 2 and m = 0 are shown. The coefficients were computed from full T-matrix (open circles), in the forward scattering approximation (squares) and the kinematic approximation (diamonds). (b) Comparison of far-field scattering amplitudes for N = 2, ka = 30, np = 1.01, kd = 2ka using full T-matrix (solid blue), in forward approximation (green dashed line) and kinematic approximation (red dashed–dotted). The forward approximation is indistinguishable from the exact solution in this case.

Fig. 7[link](b) shows a comparison of the computed far-field scattering amplitudes for N = 2, ka = 30, np = 1.01, kd = 2ka 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 = 3ka 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[link](a) and 8[link](b) show the evolution of the error [equation (14)[link]] 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.

[Figure 8]
Figure 8
(a) Error [equation (14)[link]] on the scattering amplitude coefficients with increasing number of spheres N for ka = 7, np = 1.001, kdka = 3. Squares: forward scattering approximation; diamonds: kinematic approximation; triangles: secondary approximations; coloured dashed lines: successive n-times multiple scattering approximations. (b) Error on the scattering amplitude coefficients in the increasing successive approximations for ka = 7, kdka = 3 and N = 100 spheres while varying np.

The forward scattering approximation becomes worse for increasing numbers of spheres as a result of accumulated errors. This is a limitation of a forward propagation scheme such as the one used in the MS approach, whose accuracy gets worse with increasing simulated sample thickness (Kirkland, 2019[Kirkland, E. J. (2019). Advanced Computing in Electron Microscopy, 3rd ed. Cham: Springer.]). 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)[link], (5)[link] 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[Cowley, J. M. & Moodie, A. F. (1957). Acta Cryst. 10, 609-619.]) that the scattered beam could be expressed using expansion (17)[link], considering only the scattering terms of the primary beam [n = 1,3,5\ldots]. The famous two-beam analytical scattering expression is obtained as the sum of the infinite series.

Fig. 9[link](a) shows the scattering probabilities of n-times scattering for ka = 7, np = 1.01, ka = 3kd with increasing N. The dynamic scattering probability [P_{\rm dyn} = \sum_{n\gt 1}P_{n}] 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 Poisson distribution (Egerton, 2011[Egerton, R. (2011). In Electron Energy-loss Spectroscopy in the Electron Microscope. New York: Springer.]), also shown as dashed–dotted lines calculated using the average scattering cross section. The main reason for these trends is that the ballistic picture ignores the effect of interference, which drastically affects the overall scattering process.

[Figure 9]
Figure 9
(a) Scattering probabilities (dashed coloured lines) of n-times scattering, kinematic (cyan solid line), dynamic scattering (blue solid line) and total scattering (black solid) probabilities. The total scattering probability following the Poisson distribution using the average scattering cross section is shown as the black dashed–dotted line. Probabilities have been integrated over angular distribution using ka = 7, np = 1.01 and ka = 3kd. (b) Scattering amplitudes fn, associated with n-times scattering (for n = 1-, 2- and 3-level scattering) for ka = 7, np = 1.01 and N = 15. The total scattering amplitude [f = \sum_{n}f_{n}] is shown in solid black.

Fig. 9[link](b) shows the scattering amplitudes fn, associated with n-times scattering for ka = 7, np = 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 cross section 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 transmission electron microscopy. Decreasing the energy would result in lower normalized radii, which would correspond to shifting the blue curve for carbon in Fig. 5[link](a) to the left. As a result, spherical harmonics of lower orders would be sufficient to achieve good accuracy. On the other hand, the forward scattering approximation would lose accuracy at smaller thicknesses, and then full inversion of the system (6)[link] 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 dynamical diffraction effects thanks to its ability to solve the wave equation exactly. For a 200 keV electron beam, each spherical coefficient ap, bp 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 2n2N since azimuthal coefficients [m\neq 0] 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)[link], (16)[link] 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[Kohn, W. & Rostoker, N. (1954). Phys. Rev. 94, 1111-1120.]).

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[Natoli, C. R., Benfatto, M. & Doniach, S. (1986). Phys. Rev. A, 34, 4682-4694.]). 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 structure determination 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 inelastic scattering, which can provide a greater insight into dynamical diffraction effects in real experiments. Incoherent scattering could, for example, be included in a stochastic fashion as randomly affecting the phase of the scattered waves. Inelastic scattering might require the use of dedicated electrostatic absorption potentials. It is anticipated that incoherent and inelastic scattering may have a significant effect in dynamical diffraction 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 [a_{p\semi lm}], [b_{p\semi lm}] are found by imposing the continuity of the wavefunction and its radial derivative at the surface of each sphere p:

[\left(\textstyle\sum\limits_{q = 1}^{N}\Psi_{q}^{({\rm out})}+\Psi^{({\rm i})}\right)\big{|}_{r _{p} = a_{p}} = \left(\Psi_{p}^{({\rm in})}\right)\big{|}_{r_{p} = a_{p}}]

[\partial_{r_{p}}\left(\textstyle\sum\limits_{q = 1}^{N}\Psi_{q}^{({\rm out})}+\Psi^{({\rm i})}\right) \big{|}_{r_{p} = a_{p}} = \partial_{r_{p}}\left(\Psi_{p}^{({\rm in})}\right)\big{|}_{r_{p} = a_{p}},]

where [\Psi^{({\rm i})}] is the incident electron wavefunction, [\Psi_{p}^{({\rm in})}] the scattered wavefunction of sphere p expressed inside the sphere, [\Psi_{p}^{({\rm out})}] the scattered wavefunction of sphere p expressed outside sphere p, and ap the radius of sphere p.

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

[\eqalignno{ a_{p\semi lm}& = u_{p\semi l}c_{lm}+u_{p\semi l}\textstyle\sum\limits_{q\neq p}^{N}&\cr &\quad\textstyle\sum\limits_{\nu = 0}^{\infty}\sum\limits_{\mu = -\nu}^{\mu = \nu}a_{\nu,\mu\semi l,m}^{({\rm out}{\hbox{-}}{\rm in})}({\bf d}_{pq})b_{q\semi \nu\mu}&(15)}]

[\eqalignno{ b_{p\semi lm}& = v_{p\semi l}c_{lm}+v_{p\semi l}\textstyle\sum\limits_{q\neq p}^{N}&\cr &\quad\textstyle\sum\limits_{\nu = 0}^{\infty}\sum\limits_{\mu = -\nu}^{\mu = \nu}a_{\nu,\mu\semi l,m}^{({\rm out}{\hbox{-}}{\rm in})}({\bf d}_{pq})b_{q\semi \nu\mu},&(16)}]

where the translational addition theorem (Dufva et al., 2008[Dufva, T. J., Sarvas, J. & Sten, J. C. (2008). Prog. Electromagn. Res. B, 4, 79-99.]) has been used to express the field, scattered by sphere q in the reference frame of sphere p, formally written as [\Psi_{q}^{({\rm out})}({\bf r}_{p})]. This operation involves the coupling coefficients [a_{\nu,\mu\semi l,m}^{({\rm out}{\hbox{-}}{\rm in})}({\bf d}_{pq})], where [{\bf d}_{pq} = {\bf d}_{q}-{\bf d}_{p}].

The coefficients clm are related to the incident wave. In the case of a plane wave [\exp(j{\bf k}_{0}\cdot{\bf r})], [j = \sqrt {-1} ], the addition theorem is used to expand the plane wave on the spherical Bessel solutions basis set:

[\eqalign{ c_{lm} &= 4\pi j^{l}Y_{l}^{m*}\left(\theta_{\rm i},\phi_{a}\right)e_{p} \cr e_{p}& = \exp(jk_{0}d_{p}\zeta_{p})\cr \zeta_{p} &= \sin(\Theta_{p})\sin(\Phi_{p})\sin(\theta_{\rm i})+\cos(\Theta_{p})\cos(\theta_{\rm i}),}]

where [d_{p},\Theta_{p},\Phi_{p}] are the spherical coordinates of the centre of sphere p in the global coordinate system and [0\leq\theta_{\rm i}\leq\pi] is the angle of incidence with respect to the [{\bf e}_{z}] axis. For propagation in the (y,z) plane, assumed in this paper, [\phi_{\rm i} = \pi/2] and ep is the phase offset at sphere p. The different notations are illustrated in Fig. 1[link].

The coefficients [u_{p\semi l}] and [v_{p\semi l}] are expressed as

[u_{p\semi l} = {{h_{l}^{\prime}(k_{0}a_{p})j_{l}(k_{0}a_{p})-h _{l}(k_{0}a_{p})j_{l}(k_{0}a_{p})^{\prime}} \over {j_{l}(k_{p}a_{p})h_{l}^{ \prime}(k_{0}a_{p})-n_{p}j_{l}^{\prime}(k_{p}ad_{p})h_{l}(k_{0}a_{p})}}]

[v_{p\semi l} = {{n_{p}j_{l}^{\prime}(k_{p}a_{p})j_{l}(k_{0}a_{ p})-j_{l}(k_{p}a_{p})j_{l}^{\prime}(k_{0}a_{p})} \over {j_{l}(k_{p}a_{p})h_{l}^{ \prime}(k_{0}a_{p})-n_{p}j_{l}^{\prime}(k_{p}a_{p})h_{l}(k_{0}a_{p})}},]

where [z_{l}^{\prime} = \partial_{\rho}z_{l}(\rho)].

Each of equations (15)[link], (16)[link] represent n2 equations for individual spheres where n = lmax is the truncation longitudinal order. Note that each equation involves the outside coefficients [b_{q\semi l,m}] of all the other spheres through the translational coefficients. As a result, equations (15)[link], (16)[link] represent a N×2×n2 system with as many unknowns as the number of equations and, therefore, can be written in the matrix form (6)[link].

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:

[\eqalignno{ f(h,k)&\propto\textstyle\sum\limits_{l}\sum\limits_{h_{1}}\sum\limits_{k_{1}}\sum\limits_{l_{1}}\ldots\sum\limits_ {h_{N-1}}\sum\limits_{k_{N-1}}\sum\limits_{l_{N-1}}Q_{h_{1},k_{1},l_{1}}\ldots &\cr &\quad\times Q_{h_{N-1},k_{N-1},l_{N-1}}Q\left(h-h_{u},k-k_{v},l-l_{w}\right )&\cr &\quad\times \exp\left[-2\pi j\left(H\zeta-\Delta z\textstyle\sum\limits_{n = 1}^{N- 1}\zeta_{n}\right)\right],&(17)}]

where [h = \sum_{n = 1}^{N}h_{n}], [k = \sum_{n = 1}^{N}k_{n}], [h_{u} = \sum_{n = 1}^{N-1}h_{n}], [k_{v} = \sum_{n = 1}^{N-1}k_{n}], [l_{w} = \sum_{n = 1}^{N-1}l_{n}], [H = N\Delta z] is the total thickness of the sample, N the number of slices of thickness [\Delta z], Qh,k,l = [-j/\Delta z\delta_{h,k}\exp(-2j\pi l_{n}z_{n}/c)] + [\sigma F_{h,k,l}], Fh,k,l is the structure factor, [\sigma = 2\pi m_{\rm e}h/\lambda] the interaction parameter, ζ the excitation error of beam (h,k,l) and [\zeta_{n}] is the excitation error of beam [(\sum_{r = 1}^{n}h_{r},\sum_{r = 1}^{n}l_{r})]. The excitation error is expressed as

[\zeta = {{1} \over {2K}}\left({{h} \over {a^{2}}}+{{k} \over {b^{2}}}\right)-{{l} \over {c}},\eqno(18)]

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

[Figure 10]
Figure 10
Multiple scattering in MS. (a) Distance [\zeta^{({\rm MS})}] to the Ewald paraboloid (solid curve) as represented by MS and distance [\zeta^{({\rm BW})}] to the Ewald sphere (dashed curve), known as the excitation error in the Bloch wave theory. Point P is the projection of reciprocal point (un,wn) onto the Ewald paraboloid. (b) Multiple scattering in reciprocal space for N = 3 slices located at z1,z2,z3. In this example, the beam is scattered in the first slice at z1 in direction (h1,0) = (2,0). Then, at slice z2, it scatters in direction (h2,0) = (1,0), which results in an overall scattering of the original beam in (h1+h2,0) = (3,0). Finally, at the last slice zN = z3, deflection is in direction (hN,0) = (2,0). Therefore, the overall scattering for this example corresponds to a contribution to reflection (h1+h2+h3,0) = (5,0) in the diffraction pattern. The open blue circles show how to interpret the successive excitation errors [\zeta_{i}].

From equation (17)[link] it is possible to gather terms in powers of [Q_{0,0,0}^{N-n}] into fh,k(n), which corresponds to the incident beam scattered n times before contributing to reflections h,k. The term Q0,0,0N only appears for h = k = 0 and corresponds to the unscattered direct beam. There are N terms involving the factor [Q_{0,0,0}^{N-1}] depending on which slice the single scattering event took place. There are [N(N-1)/2] terms involving the factor [Q_{0,0,0}^{N-2}] 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[Cowley, J. M. & Moodie, A. F. (1957). Acta Cryst. 10, 609-619.]). This multiple scattering picture is illustrated in 2D in Fig. 10[link](b) for a case with N = 3 and using only beams in the zero-order Laue zone (ZOLZ) li = 0.

Using only ZOLZ beams, fh,k(1), [(h,k)\neq(0,0)] is expanded as

[f_{h,k}^{(1)}\propto F_{h,k,0}\textstyle\sum\limits_{m = 1}^{N}\exp(-2j\pi m\Delta z\zeta),\eqno(19)]

which is the well known kinematic scattering regime, where the sum converges to the standard Ewald sphere curvature factor [\sin(\pi\zeta H)/\pi\zeta] for infinitely thick slices [N\rightarrow\infty], [\Delta z\rightarrow 0], [N\Delta z = H].

For a two-level scattering using only ZOLZ beams, fh(2) would expand as (written in 2D here)

[\eqalignno{ f_{h}^{(2)}&\propto\textstyle\sum\limits_{h_{1}}F_{h_{1},0}F_{h-h_{1},0}\sum\limits_{m_{1 } = 1}^{N}\sum\limits_{m_{2}\gt m_{1}}^{N}&\cr &\quad\exp\{-2j\pi\Delta z[m_{1}\zeta_{1}+m_{2} (\zeta-\zeta_{1})]\}.&(20)}]

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

First citationAllen, L. J., D'Alfonso, A. J. & Findlay, S. D. (2015). Ultramicroscopy, 151, 11–22.  Web of Science CrossRef CAS PubMed Google Scholar
First citationBalanis, C. (1989). Advanced Engineering Electromagnetics. New York: John Wiley and Sons.  Google Scholar
First citationBethe, H. A. (1928). Ann. Phys. 392, 55–129.  CrossRef Google Scholar
First citationChen, J. H., Van Dyck, D. & Op de Beeck, M. (1997). Acta Cryst. A53, 576–589.  CrossRef CAS Web of Science IUCr Journals Google Scholar
First citationChuang, 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
First citationCowley, J. M. (1995). Diffraction Physics, 3rd ed. Amsterdam: North Holland Publishing Company.  Google Scholar
First citationCowley, J. M. & Moodie, A. F. (1957). Acta Cryst. 10, 609–619.  CrossRef IUCr Journals Web of Science Google Scholar
First citationDederichs, P. H. (1971). Berichte der Kernforschungsanlage Jülich No. 797 (printed as manuscript).  Google Scholar
First citationDrevon, T. R. (2021). pyScatSpheres. https://pypi.org/project/pyScatSpheres/Google Scholar
First citationDufva, T. J., Sarvas, J. & Sten, J. C. (2008). Prog. Electromagn. Res. B, 4, 79–99.  CrossRef Google Scholar
First citationEgel, 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
First citationEgerton, R. (2011). In Electron Energy-loss Spectroscopy in the Electron Microscope. New York: Springer.  Google Scholar
First citationEremin, J. A., Orlov, N. V. & Rozenberg, V. I. (1995). J. Atmos. Terrestrial Phys. 57, 311–319.  CrossRef Web of Science Google Scholar
First citationGemmi, 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
First citationGlaeser, R. M. & Downing, K. H. (1993). Ultramicroscopy, 52, 478–486.  CrossRef CAS PubMed Web of Science Google Scholar
First citationGodin, O. A. (2011). J. Acoust. Soc. Am. 130, EL135–EL141.  Web of Science CrossRef PubMed Google Scholar
First citationHamid, A.-K., Ciric, I. R. & Hamid, M. (1990a). Can. J. Phys. 68, 1419–1428.  CrossRef Web of Science Google Scholar
First citationHamid, A.-K., Ciric, I. R. & Hamid, M. (1990b). Can. J. Phys. 68, 1157–1165.  CrossRef Web of Science Google Scholar
First citationHirsch, P. B., Howie, A., Nicholson, R. B., Pashley, D. W. & Whelan, M. J. (1965). Electron Microscopy of Thin Crystals. London: Butterworths.  Google Scholar
First citationHowie, A. (1963). Proc. R. Soc. Lond. A, 271, 268–287.  CAS Google Scholar
First citationHowie, A. (2014). J. Phys. Conf. Ser. 522, 012001.   Google Scholar
First citationIshizuka, K. (1982). Acta Cryst. A38, 773–779.  CrossRef CAS Web of Science IUCr Journals Google Scholar
First citationIshizuka, K. & Uyeda, N. (1977). Acta Cryst. A33, 740–749.  CrossRef CAS IUCr Journals Web of Science Google Scholar
First citationKirkland, E. J. (2019). Advanced Computing in Electron Microscopy, 3rd ed. Cham: Springer.  Google Scholar
First citationKlar, P. B., Xu, H., Steciuk, G., Cho, J., Zou, X., Palatinus, L. & Republic, C. (2021). chemrxiv, pp. 1–25.  Google Scholar
First citationKohn, W. & Rostoker, N. (1954). Phys. Rev. 94, 1111–1120.  CrossRef CAS Web of Science Google Scholar
First citationKorringa, J. (1947). Physica, 13, 392–400.  CrossRef Web of Science Google Scholar
First citationKorringa, J. (1994). Phys. Rep. 238, 341–360.  CrossRef Web of Science Google Scholar
First citationKottke, A. (2020). pygmm. https://zenodo.org/badge/latestdoi/21452/arkottke/pygmmGoogle Scholar
First citationLatychevskaia, T. & Abrahams, J. P. (2019). Acta Cryst. B75, 523–531.  Web of Science CrossRef IUCr Journals Google Scholar
First citationMetherell, A. & Fisher, R. (1969). Phys. Status Solidi B, 32, 551–562.  CrossRef CAS Web of Science Google Scholar
First citationMoine, O. & Stout, B. (2005). J. Opt. Soc. Am. B, 22, 1620.  Web of Science CrossRef Google Scholar
First citationMüller, H. J. (1965). Physica, 31, 688–692.  Google Scholar
First citationNannenga, B. L. & Gonen, T. (2019). Nat. Methods, 16, 369–379.   Web of Science CrossRef CAS PubMed Google Scholar
First citationNannenga, B. L., Shi, D., Leslie, A. G. W. & Gonen, T. (2014). Nat. Methods, 11, 927–930.  Web of Science CrossRef CAS PubMed Google Scholar
First citationNatoli, C. R., Benfatto, M. & Doniach, S. (1986). Phys. Rev. A, 34, 4682–4694.  CrossRef CAS PubMed Web of Science Google Scholar
First citationOleynikov, P., Hovmöller, S. & Zou, X. D. (2007). Ultramicroscopy, 107, 523–533.   Web of Science CrossRef PubMed CAS Google Scholar
First citationOphus, C. (2017). Adv. Struct. Chem. Imaging, 3, 1–11.  CrossRef PubMed Google Scholar
First citationPalatinus, L., Jacob, D., Cuvillier, P., Klementová, M., Sinkler, W. & Marks, L. D. (2013). Acta Cryst. A69, 171–188.  Web of Science CrossRef CAS IUCr Journals Google Scholar
First citationRehr, 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
First citationRehr, 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
First citationSchiff, L. I. (1955). Quantum Mechanics, 2nd ed. New York: McGraw-Hill.  Google Scholar
First citationSébilleau, D., Gunnella, R., Wu, Z. Y., Matteo, S. D. & Natoli, C. R. (2006). J. Phys. Condens. Matter, 18, R175–R230.   Google Scholar
First citationSébilleau, D., Hatada, K. & Ebert, H. (2017). Multiple Scattering Theory for Spectroscopies. Cham: Springer.  Google Scholar
First citationSébilleau, D., Natoli, C., Gavaza, G. M., Zhao, H., Da Pieve, F. & Hatada, K. (2011). Comput. Phys. Commun. 182, 2567–2579.   Google Scholar
First citationSilva, G. T., Baggio, A. L., Lopes, J. H. & Mitri, F. G. (2012). arXiv:1210.2116.  Google Scholar
First citationSubramanian, G., Basu, S., Liu, H., Zuo, J. & Spence, J. C. H. (2015). Ultramicroscopy, 148, 87–93.  Web of Science CrossRef CAS PubMed Google Scholar
First citationVainshtein, B. (1964). Structure Analysis by Electron Diffraction, pp. 114–204. Oxford: Pergamon Press.  Google Scholar
First citationVan Dyck, D. & Op de Beeck, M. (1996). Ultramicroscopy, 64, 99–107.   CrossRef CAS Web of Science Google Scholar
First citationZabloudil, J., Hammerling, R., Weinberger, P. & Szunyogh, L. (2005). Editors. Electron Scattering in Solid Matter. Berlin: Springer.  Google Scholar
First citationZuo, 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.

Journal logoFOUNDATIONS
ADVANCES
ISSN: 2053-2733
Follow Acta Cryst. A
Sign up for e-alerts
Follow Acta Cryst. on Twitter
Follow us on facebook
Sign up for RSS feeds