Potentials of individual atoms by convergent beam electron diffraction

In convergent beam electron diffraction (CBED) on two-dimensional (2D) materials, the intensity distributions within the individual CBED spots map the local atomic arrangements within the probed region. In this study we demonstrate that the average intensities within the CBED spots essentially depend on the scattering parameters by a single atom, thus, offering the possibility of the direct measurement of such parameters. Scattering potential of an individual atom can be approximated by the Gaussian function and its parameters, the standard deviation and the maximal phase shift, can be recovered from the ratio of the intensities in the zero- and higher orders spots in CBED pattern. In order to demonstrate this, we simulated CBED patterns and extracted the atomic scattering parameters from such patterns. In the simulated examples, no weak phase object approximation is applied and the proposed method provides accurate results even for materials with large phase shifts up to 1.5 rad, as for example tungsten. The scattering parameters recovered from the experimental CBED patterns of graphene and twisted graphene-hBN structure show good agreement with the theoretically obtained values.


Introduction
Quantitative estimation of electron-atom interaction is the key for understanding and interpreting electron images of biological and solid state samples. The whole analysis of biological samples in cryo-electron microscopy is built on using contrast transfer function, which is only valid in the weak phase object approximation (WPOA). The latter holds for light atoms constituting biological samples (H, C, O, N) at typical electron microscope energies (80-300 keV). The WPOA also holds for thin two-dimensional (2D) crystals such as graphene or hexagonal boron nitride (hBN). "Weak phase" means that the electron passing by atomic potential will gain an additional phase shift that is much smaller than 1 rad and the transmission function of the sample can be expanded into a Taylor series: t(x, y) = exp[iφ(x, y)] ≈ 1 + iφ(x, y), leading to linear dependency of the image contrast on the phase shift, which significantly simplifies the image analysis. The phase shift is given by φ(x,y) = σv z (x, y), where σ is the interaction parameter and v z (x, y) is the projected potential of the atom. In simulations, v z (x, y) is typically modelled using parameters tabulated in the literature [1], but little attention has been paid to recovery of the exact distribution of v z (x, y) from experimental images. Electron ptychography, a state-of-art resolution-record electron imaging technique, is currently celebrated for achieving a resolution record of 0.4 Å, that is only limited by thermal vibration of atomic lattice [2,3]. However, quantitative estimation of the phase shifts of the recovered atomic distribution has not been addressed in detail. In fact, the following controversy exists: on one hand, when theory predicts a weak phase shift, the WPOA can be applied, and a simplified data analysis is performed leading to the reconstructed phase shifts that are weak. On the other hand, the reconstructed atomic phase shifts could fall into category of weak phases only because the WPOA and the corresponding simplification were made during the reconstruction procedure. So, the question is, what are the true phase shifts of individual atoms, and, as a consequence, what are the projected potentials of such atoms? Another important questions is, whether both the phase shift and projected potential can be reconstructed directly from a diffraction pattern without a need for applying the WPOA which might bias the output results? In this study we use a 2D graphene as a test sample to investigate the electron-matter interaction. We quantitatively evaluate the phase shifts and projected potential of individual atoms from the simulated and experimental CBED pattern and compare the obtained values to the theoretical values.

CBED arrangement
Convergent beam electron diffraction (CBED) has routinely been applied for three-dimensional (3D) crystals since its invention in 1939 [4]. Recently, CBED has been demonstrated for 2D crystals, where the data analysis is straight forward and the information about lattice defects [5,6], and interlayer distances in bilayer structures can be directly retrieved from interference patterns in individual CBED spots [6,7]. A diagram for a CBED experiment is shown in Fig. 1a. In conventional electron diffraction a parallel electron beam is probing the sample and the resulting diffraction pattern exhibits sharp Bragg peaks which occupy a few pixels on the detector. The proposed here analysis requires recording of the exact intensity values of all diffraction peaks, including the zero-order peak, which is nearly impossible to realize in a conventional diffraction mode. Firstly, because the intensity values of individual sharp diffraction peaks cannot be accurately measured due to the finite pixel size, so that even the diffraction peaks of the same order would exhibit different intensity. (2) Secondly, the zero-order peak intensity is rarely accessible in an electron diffraction experiment because it is blocked due its high intensity. In CBED, a convergent electron beam is probing the sample which turns the sharp diffraction peaks into finite size spots. The other advantage of CBED of 2D materials over conventional diffraction is that the intensity distribution in a CBED spot can be directly mapped to the atomic distribution in the probed region. When CBED disks are not overlapping, the intensity distribution within a CBED disk can be analysed as a hologram and the 3D atomic positions can be retrieved. A CBED with non-overlapping disks is acquired when the convergence semiangle α is selected to be smaller than half of the Bragg diffraction angle (19.6 mrad for graphene). In addition, for quantitative analysis, presented in this paper, the detector settings should be selected in such a way, that all CBED disks are simultaneously recorded in one CBED pattern, that is, including the zero-order CBED spot.
For a perfectly flat uniform lattice, intensities within the CBED spots are constant, which we will refer to as an "average intensity in a CBED spot": where I (n) (k x , k y ) is the intensity value at (k x , k y ) pixel in (n)-th CBED spot, and A (n) is the total area of the CBED spots of (n)-th order in pixels.

Atoms as point scatterers
In the simplest model, when an electron wave illuminates the sample, each atom acts as a point scatterer and emits a spherical wave; the corresponding simulated CBED pattern is shown in Fig. 1b. The graphene lattice can be represented as sum of two trigonal lattices with a relative shift of Δx = a 0 = 142 pm, where a 0 is the closest atom-to-atom distance in graphene. The wavefront in the far field is given by: where the Fourier transform is defined as: and l 0 (x, y) is the exit wave behind one trigonal lattice. The intensity is given by where I 0 (k x , k y ) is the distribution of intensity for one trigonal lattice monolayer. From Eq.

Electron scattering on atomic potential
A more realistic model of electron scattering on a sample utilizes scattering of electrons by the potentials of individual atoms in 2D crystal. In this model, the total transmission function of the sample is written as where φ(x, y) is a phase shift by a single atom in the lattice, and l(x, y) is the function describing the positions of individual atoms in the lattice.
The sample phase shift is given by Φ(x, y) = σV z (x, y), where V z (x, y) is the sample projected potential. The interaction parameter σ is given by σ = 2πmeλ where m is the relativistic mass of the electron, e is an elementary charge, λ is wavelength of the electrons, υ is the electron velocity, and h is the Planck's constant.
where v z (x, y) is the projected potential of a single atom. Projected potential of an individual atom v z (x, y) can be modelled using parameters tabulated in the literature [1,8]. Using these, the maximal phase shift of an electron beam of 80 keV calculated at radial position of 0.1 Å amounts to: 0.217 rad for carbon and 0.214 rad for hBN (calculated as an average of 0.190 rad for B and 0.238 rad for N).

Weak phase object approximation (WPOA)
A CBED pattern of graphene simulated via scattering of electrons by the atomic potentials is shown in Fig. 1c. It was calculated using the transmission function given by Eq. (5), without approximations. One significant difference between the CBED patterns shown in Fig. 1 b and c is that the zero-order CBED spot is about 10 4 times more intense when electron scattering on the atomic potentials is taken into account. Mathematically, it can be explained as follows. For weak phase, the transmission function given by Eq. (5) can be expanded in Taylor series as where term "1" is much larger than the second term. t(x, y) = 1 corresponds to the situation without any sample and thus describes the strong signal in the zero-order CBED spot, and only in the zero-order CBED spot. The second term in Eq. (6) gives rise to a periodical signal in Fourier space, that is CBED spots of all orders, including the zero-order CBED spot; however, the signal in the zero-order CBED spot, which comes from the second term is negligibly small when comparing to the contribution from the unscattered wave described by the first term. In all simulations shown in this study, no WPOA is applied, and the transmission function is calculated by Eq. (5).

Gaussian distribution phase shift parameters evaluated from CBED pattern
The intensity within the zero-order CBED spot is given mainly by the unscattered probing wave, the contribution from the scattered wave is about 10 3 times smaller. The signals within higher-order CBED spots are given only by the scattered wave, and thus solely determined by the scattering potential of individual atoms in the lattice. Thus, neglecting the contribution from the scattered wave, the intensity of the probing beam can be determined from the averaged intensity in the zero-order CBED spot. Next, knowing the exact quantitative distribution of the probing wave, in particular its amplitude, the parameters of the scattering potential of individual atoms can be recovered from the intensities in the first-and higher order CBED spots. In the next subsections we will show that the phase shift introduced by an individual atom can be modelled as a Gaussian function and we derive the equations and demonstrate how the amplitude and standard deviation of this Gaussian function can be retrieved from the intensities of zero-and higher-orders CBED spots.

Phase shift modelled as a Gaussian distribution
Phase shift, which an electron wave acquires by interaction with an individual atom, can be accurately modelled using the parameters tabulated in Refs. [1,8]. In this study, the phase shift is modelled in the form of a Gaussian function: with amplitude φ 0 and standard deviation σ p . This modelling allows us to freely change the parameters and see whether the same values of the parameters are reconstructed by the proposed method. Fig. 1d shows the phase shift of a single carbon atom calculated using Kirkland's parameters, and a superimposed fit by a Gaussian function such as Eq. (7) with φ 0 = 0.217 rad and σ p = 0.252 Å.
For larger φ 0 , and the coordinates (x, y) in the vicinity of the atom, the following Taylor series expansion can be applied giving: indeed gives a sharp peak in the center of the spectrum due to the constant terms and it gives a Gaussiandistributed signal due to the Fourier transform of iφ 0 exp The Fourier transform of a Gaussian function is also a Gaussian function, and the Fourier transform of φ(x, y) is given by:

Finding the standard deviation
Using the Taylor series expansion (in the WPOA or as described above for strong phase shifts), the intensity distribution of a CBED pattern can be written as: where p is the complex-valued distribution of the probing wave, and FT means Fourier transform. The average intensity of the zero-order CBED spot is mainly given by the unscattered wave (term "1" in Eq. (10)): The signal in the higher-order CBED spots is provided by the scattered wave only. For a Gaussian-distributed phase shift of individual atoms, the average intensity in a higher-order CBED spot can be found via analytical solution: where k (n) is the coordinate of the center of the n-th order CBED spot and p 0 is the amplitude of the probing wave. From Eq. (12) it can be seen that the standard deviation σ p,F can be directly recovered by comparing the averaged intensities in the first-and the second order CBED spots: where the factor 4 accounts for the difference in the intensity due to the hexagonal lattice, without projected potential, as discussed above.

Finding the amplitude
The Fourier transform of the phase shift function φ(x, y) in Cartesian coordinates is given by: where The averaged intensity in the zero-th order CBED spot is given approximately by the probing wave alone and amounts to: where M = z/Δf is the magnification factor, Δf is the defocus distance, z is the sample-to-detector distance (camera length), and I 0 is the intensity of the probing wave. The averaged intensity in a first-order CBED spot of a monolayer depends on the number of the scattering atoms (that is, the number of atoms in the probed region) and on the phase shift distribution parameters as follows: where σ d is given by Eq. (15), is the coordinate of the center of the CBED spot, N 0 is the number of atoms in the probed area, 2 , h 0 is the area of one 2D unit cell (parallelogram) in the graphene lattice, and the factor 2 accounts for two carbon atoms per unit cell. Note that 〈I (0) 〉 and 〈I (1) 〉 both contain factor I 0 , which is eliminated when calculating the ratio 〈I (1) 〉/〈I (0) 〉: where γ is the number of layers. The ratio 〈I (1) 〉/〈I (0) 〉 does not depend on the intensity of the probing beam I 0 , defocus distance Δf or the sampleto-detector distance z. By using Eqs. (9) and (13), φ 0 can be determined from the known parameters α, λ, h 0 , ratio of measured intensities 〈I (1) 〉/ 〈I (0) 〉, and σ p determined from the ratio 〈I (1) 〉/〈I (2) 〉. When the parameters σ p and φ 0 are known, the number of layers γ can be evaluated by comparing the averaged intensities in the zero-and first order CBED spots using Eq. (18).

Parameters reconstructed from simulated CBED patterns
A series of CBED patterns for atomic lattices with different parameters for the phase shift of individual atoms were simulated (as described in Appendix B) and the parameters were recovered by using the equations above. CBED patterns were simulated at defocus values Δf = 0.5, 1.0, and 2.0 μm, and it was also verified that the recovered standard deviation σ p and amplitude φ 0 did not depend on the defocus value Δf.
The scattering parameters were evaluated from the simulated CBED patterns. σ p was evaluated from the ratio 〈I (1) 〉/〈I (2) 〉 using Eq. (13), and φ 0 from the ratio 〈I (1) 〉/〈I (0) 〉 using Eq. (18). The standard deviation σ p was correctly reconstructed as 0.2 Å from all the CBED patterns. The reconstructed phase shifting amplitudes φ 0 slightly deviated from the original values, with the mismatch increasing towards higher φ 0 , reaching a mismatch of 0.04 rad for φ 0 = 1.5 rad, as illustrated in Fig. 1e. This result is quite interesting. φ 0 = 1.5 rad is a large phase shift, for which the conventional WPOA cannot be applied. The corresponding CBED pattern was simulated without any approximations and the recovered (relatively large) phase shift φ 0 is almost correct. This confirms that the approximation for strong phase shifts which we derived above can be applied and leads to correct results. The proposed analysis can be applied for determination of scattering potentials of arbitrary atoms which cause even larger phase shifts than required by the conventional WPOA. This could include higher atomic number atoms, such as molybdenum (φ 0 = 0.969 rad at 80 keV) or tungsten (φ 0 = 1.415 rad at 80 keV), or for other atoms in transition metal dichalcogenide (TMD) structures.
The CBED pattern of a bilayer 2D material structure, where atoms in each layer exhibit σ p = 0.2 Å and φ 0 = 0.1 rad, with a twist angle of 2 • was simulated. In the analysis, only the areas where the CBED spots from individual layers overlap were considered (regions with overlapping and non-overlapping CBED spots are illustrated in Fig. 2b). For a bilayer system, the average intensity in a CBED spot amounts to the sum of the intensities from individual layers. σ p was determined from the ratio 〈I (1) 〉 /〈I (2) 〉 using Eqs. (9) and (13); here the intensity in both the first-and second-order CBED spots was higher by factor of 2 due to the presence of two layers. φ 0 was determined from the ratio 〈I (1) 〉/〈I (0) 〉 using Eq. (18); here the two layers led to the intensity in the first-order CBED spots increasing by factor of 2, so that the measured intensity 〈I (1) 〉 needed to be reduced by factor 2 in Eq. (18). The recovered values are σ p = 0.212 Å and φ 0 = 0.096 rad. When comparing to the original values, the errors are small, and can be explained by the intensity values being averaged over the interference patterns in the first-and second-order CBED spots, where the amount of maxima and minima might be not exactly the same so as to get the correct average.

Monolayer graphene
The experimental CBED pattern of monolayer graphene is shown in Fig. 2a. The samples were prepared and imaged in CBED mode as described in Appendix C, and the obtained CBED patterns were analysed as described in Appendix D. The following parameters of the Gaussian distribution were reconstructed by the proposed method: σ p = 0.308 Å and φ 0 = 0.156 rad, which compares well to the theoretically calculated values mentioned above: σ p = 0.252 Å and φ 0 = 0.217 rad. The main source of the mismatch is the noise, and the presence of the ripples or adsorbates on the graphene surface which lead to the intensity variations within CBED spots, as it can be seen in the CBED pattern shown in Fig. 2a.

Bi-layer graphene-hBN
The experimental CBED pattern of twisted bilayer graphene-hBN is shown in Fig. 2b. In this CBED pattern, the contribution from individual layers can be extracted from the regions of non-overlapping disks. The signal in the CBED disks can be easily assigned to the individual layers. The CBED disks observed at higher scattering angles are from the graphene layer (smaller interatomic distance), and the CBED disks observed at lower scattering angles are from the hBN (larger interatomic distance). By using the quantitative analysis, we obtain σ p = 0.229 Å and φ 0 = 0.189 rad for the graphene layer and σ p = 0.244 Å and φ 0 = 0.162 rad for the hBN layer. As predicted by theory, the maximal phase shift is Fig. 2. Experimental convergent beam electron diffraction (CBED) patterns of (a) monolayer graphene and (b) monolayer graphene on monolayer hBN. Both CBED patterns acquired at 80 keV and Δf = 1.8 μm. The zero-order spot intensity is reduced by factor 10 − 3 . In (b), a region with overlapping spots is indicated by the orange arrow, regions with non-overlapping spots are indicated by the yellow arrows. Scalebars are 2 1/nm. higher for graphene. Quantitatively, the maximum phase shifts for both layers are somewhat smaller than expected from the theory, φ 0 = 0.217 rad for graphene and φ 0 = 0.214 rad for hBN.
For a twisted bilayer system, interference pattern is observed in the regions where the CBED disks overlap [6]. The averaged intensity of such an interference pattern is given by the sum of the averaged intensities from the two layers. By using quantitative analysis on the areas where the CBED disks overlap, the following parameters were recovered σ p = 0.295 Å and φ 0 = 0.153 rad. The reconstructed standard deviation σ p appears larger than determined from the CBED intensity of individual layers and the maximal phase shift φ 0 is smaller than that recovered for individual layers. This mismatch is similar to the mismatch observed for the simulated CBED pattern of the bilayer system. Similar to the case of the simulated CBED pattern for the bilayer system, the errors here can be explained by errors introduced by averaging over the interference pattern.

Conclusions
We derived equations for recovery of single atom scattering parameters, the maximum phase shift and the extent of the scattering potential, directly from corresponding CBED pattern and ratio of the intensities in the zero-and higher orders CBED spots. The scattering parameters recovered from simulated and experimental CBED patterns exhibit good agreement with the theoretical values. The simulations confirm that the proposed analysis can also be applied for atoms which scatter strongly and introduce relatively large phase shifts into the probing electron beam. Thus, the proposed data analysis could be utilized for determination of scattering potentials of arbitrary atoms, including atoms in TMD structures.
The main sources of errors are the noise and presence of ripples or adsorbates on the graphene surface which lead to intensity variations within CBED spots. For bilayer twisted systems, the defocus values should be selected to be relatively large so that the interference patterns in the individual CBED spots exhibit several periods and the averaged intensity can be extracted accurately.
Extraction of the number of layers from the electron diffraction data proved to be a challenging task in electron diffraction imaging of van der Waals structures [9,10]. When the parameters of individual atom phase shift are known, the number of layers can be determined from the ratio of the average intensities in the zero-and first orders.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability
Data will be made available on request.

A. Fourier transform of function with Gaussian-distributed phase
The following approximation exp[iφ(x, y)] = exp
Probing wave distribution. The probing wave distribution ψ 0 ( r → ) was modelled by simulating the diffraction of a spherical wavefront on a limiting virtual aperture (second condenser aperture) positioned at a plane r → 0 : where a( r → 0 ) is the aperture binary function. In the simulations, the integral transform in Eq. (A.1) was calculated directly, without applying a fast Fourier transform (FFT), as where m, n, p, q are the pixel indices in the aperture (m, n) and the sample (p, q) planes, and Δ a and Δ are the pixel sizes in the aperture and the sample planes, respectively. The convergence semiangle of α = 7 mrad was achieved by selecting a virtual exit aperture of diameter 640 nm at a distance of 45 μm from the virtual source plane.

CBED patterns of monolayer samples
The input data consisted of an array of the coordinates of all the atoms (x n ,y n ). The transmission function of a monolayer was calculated as: where v z (x, y) is the projected potential of an individual atom, l(x, y) is the function describing positions of the atoms in the lattice, and ⊗ denotes convolution. The projected potential of a single carbon atom was simulated in the form ref. [1]: where r = ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ ̅ x 2 + y 2 √ , a 0 is the Bohr radius, e is the elementary charge, K 0 (...) is the modified Bessel function, and a i , b i , c i , d i are parameters that depend on the chemical origin of the atoms and are tabulated elsewhere [1]. In v z (r), the singularity at r = 0 is replaced by the value of v z (r) at r = 0. The exit wave after passing through the sample is given by the product of the incident wave ψ 0 (x, y) and the transmission function of the sample u(x, y) = ψ 0 (x, y)t(x, y). The CBED pattern is then simulated as the square of the amplitude of the FT of the exit wave. For a bilayer system, the transmission function of the bilayer was calculated as a product of transmission functions of individual layers.

C. Sample preparation and acquisition
The samples were obtained by mechanical stacking of mechanically exfoliated graphene and hBN layers on a Si/SiO 2 substrate. By using "pick and lift" technique [12] -the whole stack was then transferred on a quantifoil carbon support film for observation in the TEM. CBED was realized in a probe side aberration corrected scanning TEM at an accellerating voltage of 80 keV and a convergence semiangle, α, of ~6-8 mrad. During experiment, the focal lengths of the objective and condenser lenses were kept constant, thus there was no change in the convergence angle. The sample z-position was changed by moving the sample along the optical axis. Bilayer structure significantly improves the stability of 2D crystals upon exposure to high-energy electrons and the electron dose required for CBED imaging is low, so no evidence of knock-on damage was observed during prolonged data acquisition [13,14].

D. Analysis of experimental CBED patterns
For experimental images, it was observed that the noise level is decreasing towards higher scattering angles (higher k-values). To minimize the effect of the noise dependence on the scattering angle position, the signal and the noise were subtracted at a selected scattering angle, corresponding to the maximal extend of a CBED spot.