Bohmian trajectory-bloch wave approach to dynamical simulation of electron diffraction in crystal

In the Bohmian trajectory theory the trajectories are guided by the wave function. In order to obtain the wave function, we start with solving the Schrödinger equation for a periodic crystal potential. The wave function of incident electrons inside the periodic potential of a perfect crystal is expressed as a superposition of Bloch waves,

$$\begin{eqnarray}&&\psi ({\bf{r}})=\displaystyle \sum _{j}{A}^{(j)}b({{\bf{k}}}^{(j)},\,{\bf{r}})=\displaystyle \sum _{j}{A}^{(j)}\displaystyle \sum _{{\bf{g}}}{C}_{{\bf{g}}}^{(j)}\exp [{\rm{i}}({{\bf{k}}}^{(j)}+{\bf{g}})\cdot {\bf{r}}],\end{eqnarray} \tag{ 1 }$$

where ${A}^{(j)}$ and ${{\bf{k}}}^{(j)}$ are respectively the amplitude and the wave vector of the Bloch wave $b({{\bf{k}}}^{(j)},\,{\bf{r}}).$ In the time-independent Schrödinger equation,

$$\begin{eqnarray}&&\left[-\displaystyle \frac{{\hslash }^{2}}{2m}{{\rm{\nabla }}}^{2}+U({\bf{r}})\right]\psi ({\bf{r}})=E\psi ({\bf{r}}),\end{eqnarray} \tag{ 2 }$$

where $U({\bf{r}})$ is a periodic potential and can be expanded as a Fourier series,

$$\begin{eqnarray}&&U({\bf{r}})=\displaystyle \sum _{{\bf{g}}}{U}_{{\bf{g}}}\exp ({\rm{i}}{\bf{g}}\cdot {\bf{r}}),\end{eqnarray} \tag{ 3 }$$

by substituting equations (1) and (3) into equation (2), the renowned dynamical equation of Bethe [34] is obtained as

$$\begin{eqnarray}&&[{{\bf{K}}}^{2}-{({{\bf{k}}}^{(j)}+{\bf{g}})}^{2}]{C}_{{\bf{g}}}^{(j)}+\displaystyle \sum _{{\bf{h}}\ne {\bf{g}}}{U}_{{\bf{g}}-{\bf{h}}}{C}_{{\bf{h}}}^{(j)}=0,\end{eqnarray} \tag{ 4 }$$

where ${{\bf{K}}}^{2}={{\bf{k}}}_{0}^{2}+{U}_{0}$ is the electron wave vector inside the crystal and ${{\bf{k}}}_{0}$ is the electron wave vector in the vacuum. Then the vector ${{\bf{k}}}^{(j)}$ is decomposed as,

$$\begin{eqnarray}&&{{\bf{k}}}^{(j)}={\bf{K}}+{\lambda }^{(j)}\hat{{\bf{n}}},\end{eqnarray} \tag{ 5 }$$

where $\hat{{\bf{n}}}$ is the unit normal vector of the surface. By substituting equation (5) into equation (4), the equation is converted into an eigenvalue equation [35] in the form of

$$\begin{eqnarray}&&[{({\lambda }^{(j)})}^{2}{\bf{I}}+{\lambda }^{(j)}{\bf{H}}-{\bf{D}}]{{\bf{C}}}^{(j)}={\bf{0}},\end{eqnarray} \tag{ 6 }$$

where ${\bf{I}}$ is a unit matrix, ${H}_{{\bf{g}}{\bf{g}}}=2\hat{{\bf{n}}}\cdot ({\bf{K}}+{\bf{g}}),$ ${H}_{{\bf{g}}{\bf{h}}}=0,$ ${D}_{{\bf{g}}{\bf{g}}}=-{{\bf{g}}}^{2}-2{\bf{K}}\cdot {\bf{g}},$ ${D}_{{\bf{g}}{\bf{h}}}={U}_{{\bf{g}}-{\bf{h}}},$ and ${{\bf{C}}}^{(j)}$ is a column vector containing the Fourier coefficients ${C}_{{\bf{g}}}^{(j)}.$ Equation (6) can be solved by a standard method [36]. However, for high-energy electrons, it is a good approximation [36, 37] to keep only the linear term and ignore the high-order term; then the eigenvalue equation is rewritten as,

$$\begin{eqnarray}&&{\bf{D}}{{\bf{C}}}^{(j)}={\lambda }^{(j)}{\bf{H}}{{\bf{C}}}^{(j)}.\end{eqnarray} \tag{ 7 }$$

Define a matrix ${\bf{B}}={{\bf{H}}}^{-1}{\bf{D}}$ so that

$$\begin{eqnarray}&&{\bf{B}}{{\bf{C}}}^{(j)}={\lambda }^{(j)}{{\bf{C}}}^{(j)},\end{eqnarray} \tag{ 8 }$$

where ${B}_{{\bf{g}}{\bf{g}}}=-({{\bf{g}}}^{2}+2{\bf{K}}\cdot {\bf{g}})/[2\hat{{\bf{n}}}\cdot ({\bf{K}}+{\bf{g}})]$ and ${B}_{{\bf{g}}{\bf{h}}}={U}_{{\bf{g}}-{\bf{h}}}/[2\hat{{\bf{n}}}\cdot ({\bf{K}}+{\bf{g}})].$ Meanwhile, the amplitudes of Bloch waves ${A}^{(j)}$ can be obtained from

$$\begin{eqnarray}&&{A}^{(j)}={C}_{j}^{-1},\end{eqnarray} \tag{ 9 }$$

where ${C}_{j}^{-1}$ is the jth element in the first column of matrix ${{\bf{C}}}^{-1},$ which is demonstrated [38] by Sheinin et al With the coefficients solved above, the total wave function is obtained as:

$$\begin{eqnarray}&&\psi ({\bf{r}})=\displaystyle \sum _{j}{A}^{(j)}\displaystyle \sum _{{\bf{g}}}{C}_{{\bf{g}}}^{(j)}\exp [{\rm{i}}({\bf{K}}+{\lambda }^{(j)}\hat{{\bf{n}}}+{\bf{g}})\cdot {\bf{r}}].\end{eqnarray} \tag{ 10 }$$

Inelastic scattering can be also considered when extending the potential to be complex, U(r) = U^Re(r) + iU^Im(r). Thus, the coefficients ${A}^{(j)},$ ${C}_{{\bf{g}}}^{(j)}$ and the eigenvalues ${\lambda }^{(j)}$ all become complex. To calculate the inelastic scattering processes, thermal diffuse scattering [39, 40] is evaluated using the Debye–Waller factor [41–43], and other excitations [44] are accounted for by inelastic mean free path (IMFP) [45]. The imaginary part of the complex potential can be calculated by the FSCATT subroutine [46]. The wave function obtained is then used to calculate the Bohmian trajectory. In this paper, we consider two methods for the calculation of the Bohmian trajectory. The first one is based on the traditional Bohmian mechanics and the second one is derived from the ME approach.

In Bohmian mechanics, there exist both quantum waves and point-like particles guided by the wave functions. For large population of particles, the Bohmian mechanics gives exactly the same result as Copenhagen interpretation. This paper discusses the trajectory representation of electron diffraction in a crystal together with the particle velocity field.

The wave function $\psi ({\bf{r}},\,t)$ is complex so that it can be expressed in the polar form,

$$\begin{eqnarray}&&\psi ({\bf{r}},\,t)=R({\bf{r}},\,t)\exp [{\rm{i}}S({\bf{r}},\,t)/\hslash ],\end{eqnarray} \tag{ 11 }$$

where $R({\bf{r}},t)$ is the amplitude and $S({\bf{r}},t)$ is the phase. Then the time-dependent Schrödinger equation,

$$\begin{eqnarray}&&{\rm{i}}\hslash \displaystyle \frac{\partial \psi ({\bf{r}},\,t)}{\partial t}=-\displaystyle \frac{{\hslash }^{2}}{2m}{{\rm{\nabla }}}^{2}\psi ({\bf{r}},\,t)+U({\bf{r}},\,t)\psi ({\bf{r}},\,t)\end{eqnarray} \tag{ 12 }$$

is separated into the real part and the imaginary part [4] as

$$\begin{eqnarray}&&-\displaystyle \frac{\partial S}{\partial t}=\displaystyle \frac{{({\rm{\nabla }}S)}^{2}}{2m}+U-\displaystyle \frac{{\hslash }^{2}}{2m}\displaystyle \frac{{{\rm{\nabla }}}^{2}R}{R};\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {R}^{2}}{\partial t}+{\rm{\nabla }}\cdot \left({R}^{2}\displaystyle \frac{{\rm{\nabla }}S}{m}\right)=0.\end{eqnarray} \tag{ 14 }$$

The real part, equation (13), is analogous to the classical Hamilton–Jacobi equation, where ${({\rm{\nabla }}S)}^{2}/2m$ and $U$ respectively stand for kinetic energy and classical potential, and $Q({\bf{r}},t)=-{\hslash }^{2}{{\rm{\nabla }}}^{2}R/2mR$ is interpreted as quantum potential, which is nonlocal. The particle velocity can thus be defined as,

$$\begin{eqnarray}&&{\bf{v}}=\displaystyle \frac{{\rm{\nabla }}S}{m}=\displaystyle \frac{\hslash }{m}\text{Im}\displaystyle \frac{{\rm{\nabla }}\psi }{\psi }.\end{eqnarray} \tag{ 15 }$$

Notice that ${R}^{2}$ represents the probability density $\rho ({\bf{r}},\,t)={| \psi ({\bf{r}},t)| }^{2}={R}^{2},$ therefore, the equation (14) for the imaginary part can be rewritten as the continuity equation,

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}+{\rm{\nabla }}\cdot (\rho {\bf{v}})=0,\end{eqnarray} \tag{ 16 }$$

where ${\bf{J}}=\rho {\bf{v}}$ is the probability flow. This equation was firstly introduced in the hydrodynamic interpretation of quantum mechanics by Madelung [47]. The continuity equation here guarantees that the information of the quantum system can be computed from quantum trajectories at any time.

When simulating electron diffraction in a crystal, the wave function is obtained as equation (10), thus with the definition in equation (15) the velocity field of electrons can be calculated as,

$$\begin{eqnarray}&&{\bf{v}}=\displaystyle \frac{\hslash }{m}\text{Im}\left(\displaystyle \frac{{\sum }_{j}\displaystyle {A}^{(j)}{\sum }_{{\bf{g}}}\displaystyle {C}_{{\bf{g}}}^{(j)}{\rm{i}}({\bf{K}}+{\lambda }^{(j)}\hat{{\bf{n}}}+{\bf{g}})\exp [{\rm{i}}({\bf{K}}+{\lambda }^{(j)}\hat{{\bf{n}}}+{\bf{g}})\cdot {\bf{r}}]}{\displaystyle {\sum }_{j}{A}^{(j)}{\sum }_{{\bf{g}}}\displaystyle {C}_{{\bf{g}}}^{(j)}\exp [{\rm{i}}({\bf{K}}+{\lambda }^{(j)}\hat{{\bf{n}}}+{\bf{g}})\cdot {\bf{r}}]}\right).\end{eqnarray} \tag{ 17 }$$

Then the Bohmian trajectories are calculated by integral of the velocity.

In Bohmian mechanics, the amplitude information and the phase information of the total wave function are all contained in the Bohmian trajectories. This means that the trajectories obtained should be able to describe the details that how electrons are dynamically diffracted in a crystalline material. The Bohmian velocity calculated in periodic potential is a function of the particle position r and the wave vector K inside the crystal, which indicates that electrons with different incident directions and energies follow different trajectories. Besides the simple incident cases, this method can also deal with some complex situations when moving electrons are at different incident conditions. As an example, the simulation of Bohmian trajectories in EBSD will be further discussed in section 3.2.

In EBSD simulation, what we are mainly concerned about is how the trajectories distribute on the phosphor screen rather than how they propagate from the crystal to the screen. Accordingly, we have developed another approach named ME approach. Without performing the velocity integration, this approach gives the same trajectory distribution as the traditional Bohmian mechanics method introduced above. Therefore, it is better in efficiency and computation cost.

In the ME approach the MEs of different diffracted beams are calculated, in this way the final trajectory distribution can be directly obtained. The wave function is transformed into momentum space,

$$\begin{eqnarray}&&\varphi ({\bf{p}})=\displaystyle \frac{1}{{(2\pi \hslash )}^{3/2}}\displaystyle \int {\rm{d}}{\bf{r}}\exp (-{\rm{i}}{\bf{k}}\cdot {\bf{r}})\psi ({\bf{r}}).\end{eqnarray} \tag{ 18 }$$

Based on the wave function in the form of equation (10) we have,

$$\begin{eqnarray}&&\varphi ({\bf{p}})={(2\pi \hslash )}^{3/2}\displaystyle \sum _{j}{A}^{(j)}\displaystyle \sum _{{\bf{g}}}{C}_{{\bf{g}}}^{(j)}\delta ({{\bf{k}}}^{(j)}+{\bf{g}}-{\bf{k}});\end{eqnarray} \tag{ 19 }$$

then the ME is accordingly obtained from the integral,

$$\begin{eqnarray}&&\begin{array}{l}\langle {\bf{p}}\rangle =\displaystyle \int {\varphi }^{\ast }({\bf{p}}){\bf{p}}\varphi ({\bf{p}}){\rm{d}}{\bf{p}}\\ \,\,\,\,\,=\displaystyle \frac{1}{{(2\pi \hslash )}^{3}}\displaystyle \sum _{i,j}{A}^{(i)* }{A}^{(j)}\displaystyle \sum _{{\bf{g}},{\bf{h}}}{C}_{{\bf{g}}}^{(i)* }{C}_{{\bf{h}}}^{(j)}\delta ({{\bf{k}}}^{(i)}+{\bf{g}}-{{\bf{k}}}^{(j)}-{\bf{h}})({{\bf{k}}}^{(i)}+{\bf{g}}).\end{array}\end{eqnarray} \tag{ 20 }$$

The averaged velocity is defined as $\langle {\bf{v}}\rangle =\langle {\bf{p}}\rangle /m,$ which describes the average electron diffraction process in the lattice. With these equations the trajectory distribution in EBSD can be evaluated.

For the sake of observation and analysis, we first consider the case of normal incidence along the negative x-axis direction in the crystal (figure 1). Ideally, the Bloch wave method can yield accurate results if many reflected beams are included. In practice, a certain number of reflected beams should be selected to avoid very large matrices, which consume calculation time and may cause numerical instabilities, based on the criteria [49, 50] using the Bethe perturbation scheme. LAPACK is utilized here to solve the eigenvalue equation, equation (8), and to obtain the wave function.

Zoom In Zoom Out Reset image size

After all the coefficients in equation (10) are solved, the velocity field is acquired from equation (17). Figure 2 shows how the velocity field evolves along electron trajectory. The velocity mainly follows the electron incident direction. Figure 2(b) displays specifically the velocity field in the x–y plane at the bottom face of the lattice (i.e. the green rectangle area in figure 2(a)). The component of the velocity along the incident direction, ${{\bf{v}}}_{\parallel },$ is much larger than the perpendicular component, ${{\bf{v}}}_{\perp },$ because of the high energy of the incident electrons.

Zoom In Zoom Out Reset image size

To understand how the velocity field evolves over time, we now neglect ${{\bf{v}}}_{\parallel }$ and focus on the distribution of vertical velocity, ${{\bf{v}}}_{\perp },$ in y–z plane at different x-positions. Considering the translational symmetry, the vertical velocity distribution can be described in a single unit cell. In figures 3(a)–(d), four y–z planes (red squares) are selected in a cell (the 11st cell counted from the initial incident unit cell). On these planes the normalized vertical velocity field is shown by arrows, and the contour map represents the probability density calculated from the wave function; the velocity is seen regularly distributed under the modulation of the lattice potential. In principle, the continuity equation, equation (16), governs the relationship between the probability density ρ and the velocity field v. In figure 3, such relationship is expressed as the regular and symmetrical distribution of ρ and ${{\bf{v}}}_{\perp }.$ Simpler wave functions are calculated by reducing the number of the Fourier expansion items in figures 3(e) and (f), where this relationship is easier to observe.

Zoom In Zoom Out Reset image size

To further explore how parallel electron beams are diffracted in a crystal, the Bohmian trajectories are simulated. Assuming the uniform distribution of the initial electron trajectories landing on the crystal surface, the Bohmian trajectories are calculated from the integral of the velocity field. Figure 4(a) is a schematic of our model which consists of 50 unit cells. The blue arrows indicate that the main propagation direction of the Bohmian trajectories is along the negative x-axis inside a crystal while the positive x-axis represents the vacuum region. The trajectories are diffracted into different directions and will form a regular distribution in any y–z section when passing through the lattice, as is shown in figure 4(b). The distributions in y–z planes illustrate the proportional relationship between the probability density and the trajectory distribution. The evolution of the trajectories in 3D space is displayed in figure 4(d). It is evident that under the modulation of the lattice potential, many channels are formed along the atom column inside the crystal through which the trajectories are easier to pass. Specifically, the trajectories converge towards the channels but never intersect. In figure 4(c), which provides a side view along the incident direction, two regions are separated by a red square. Inside the square, the trajectories converge towards the center channel and the distance from the center determines the speed of the convergence, which is also indicated in figure 4(e). Similarly, outside the square the trajectories converge towards the corner channels. Different colored balls represent the sequential intersections of the trajectories with a series of y–z planes in their motion history at depth from 10 to 50 a₀, which can be easily observed from the top view displayed in figure 4 (e). Overall, the trajectories closer to the channels converge in the first few lattices, and then begin to diverge near the tenth lattice and finally converge again. While the trajectories distant from the atom columns converge more slowly. A 3D perspective in figure 4(f) illustrates the evolving of the trajectories around a channel. This motion of parallel electron beams is in fact a weakened electronic Talbot effect. The Talbot effect is that, when a plane wave is incident upon a periodic diffraction grating, the image of the grating is repeated at regular distances away from the grating plane. The Talbot effect of photons and atoms is widely known. A simulation of the atomic Talbot effect based on Bohmian mechanics is discussed in [13]. However, the electronic Talbot effect [51] is not easy to be observed in crystal because inelastic processes will weaken such an effect. Without electron inelastic scattering the trajectories will periodically repeat between convergence and divergence along the propagation direction, i.e. the Talbot effect, as is shown in figure 4(g). But under the influence of inelastic processes the electrons will finally converge and pass through the channels.

Zoom In Zoom Out Reset image size

In this simulation, we present some results regarding the diffraction of normal incident parallel electron beams calculated with the Bohmian theory by the Bloch wave method. It shows the relationship between the probability density and the velocity field as well as the trajectories. The detailed process of electron diffraction is reproduced intuitively. In particular, the weakened electronic Talbot effect and the channels of electron beams in crystal are very interesting. In theory, as an alternative formalism of quantum mechanics, the Bohmian mechanics gives a better understanding of the physical process beyond the standard formalism in which the importance of probability density is highly emphasized while phase is despised sometimes. While in this method the velocity field is directly calculated from phase, thus it presents unique advantages especially when employed to quantum coherence, interference and diffraction. Accordingly the Bohmian type of analysis can give more intuitive information by allowing us to observe the evolution of each individual trajectory with no need of a probabilistic description.

The Bloch wave-Bohmian approach will be useful to treat more complex electron diffraction when electron beams are emitted into different directions. Here we use it to simulate EBSD as an example. The existing simulation models for EBSD is summarized in [21, 52]. A dynamical simulation based on the Bloch wave theory proposed by Winkelmann et al [26] is used to compute the EBSD pattern. In the dynamical method, the intensity distribution of the diffraction pattern is obtained from the probability density. Figure 5(a) shows the simulated pattern from Mo (001) at 20 kV; the Kikuchi patterns and the higher order Laue zone (HOLZ) rings are marked by the arrows. An experimental EBSD pattern [53] is illustrated in figure 5(b). In figure 5(c), the simulated EBSD pattern and the experimental pattern are shown to agree with each other; they both serve as comparison against our trajectory approaches.

Zoom In Zoom Out Reset image size

As discussed in [54], it is a good approximation to consider that the coherence shall be destroyed by the inelastic process and, thus, the incoming and outgoing diffraction processes do not interfere. EBSD is the outgoing beam diffraction formed by electrons which are emitted into all directions from a point source inside the crystal. This process is previously simulated based on the reciprocity principle [54, 55], which states that the detected intensity of a wave on the phosphor screen originating from a point source inside the crystal, is equal to the observed intensity at the point detector inside the crystal, after diffraction of an incoming wave from the screen. Thus the outgoing beam diffraction is transformed into incoming beam diffraction. But in the present approaches we can simply and directly realize it by launching initial trajectories homogeneously into spherical directions.

The calculated Bohmian trajectories in EBSD are shown in figure 6. The electron trajectories at different emission directions are diffracted through the crystal lattice; they form a point distribution when they finally reach the screen. This is similar to the two-slit single-electron experiment [56], where an interference pattern is seen gradually formed with growing population of electrons hitting a phosphor screen. In our simulation of EBSD, by increasing the number of the trajectories gradually (here we calculate 10², 10³, 10⁴ and 10⁵ trajectories), the diffraction pattern is similarly formed from a few points on the screen to a full diffraction pattern (figures 6(a)–(e)).

Zoom In Zoom Out Reset image size

As discussed previously, there is a proportional relationship between the probability density and the trajectory distribution. This indicates that the characteristics of a diffraction pattern calculated from probability density should also be found in a trajectory distribution pattern. To find such detailed features, denser trajectories are calculated. In figure 6(e) the Kikuchi patterns and HOLZ rings appear with 10⁵ trajectories. To compare with the intensity of the simulated EBSD pattern, in figure 6(f) the black and white color is reversed to give an anti-color monochrome image which is further mixed with 30% intensity of the ESBD pattern in figure 5(a). It is shown that the trajectories are consistent with the contrast of the EBSD pattern, especially the structure of Kikuchi bands and HOLZ rings. However, this simple point representation of density distribution cannot truly display the intensity; therefore, the number of electron trajectories in each pixel of an image is counted to acquire the true statistical trajectory density distribution, which is displayed in figure 7(d).

Zoom In Zoom Out Reset image size

The results above are calculated using the Bohmian mechanics approach. We have also calculated the distribution using the ME approach. Figure 7(a) shows the point distribution of 10⁵ trajectories simulated using the ME approach; the anti-color image is displayed in figure 7(b) which presents very similar distribution with that calculated by Bohmian mechanics in figure 6(f). Then 10⁷ trajectories are calculated in order to get the trajectory density distribution or EBSD pattern using both methods: figure 7(c) using the ME approach and figure 7(d) using the Bohmian mechanics. The Kikuchi pattern and the HOLZ rings show similar features in both figures.

To measure the similarity between figures 5(a) and 7(c) and (d), the 2D image cross-correlation is calculated and shown in figure 8, where the EBSD pattern simulated by the dynamical method (figure 5(a)) is used as a reference. The image cross-correlation returns a matrix of correlation coefficients. The maximum value indicates the best relative location between the two graphs for comparison, where the two patterns most resemble each other. In our surface map of the cross-correlation matrix in figure 8, two sharp peaks appear when the dynamically simulated EBSD pattern coincides with the trajectory density distribution patterns, which quantitatively verifies that both the ME approach and the Bohmian mechanics method match well with the dynamical method. Although both the patterns of trajectory density distribution look the same, the quantitative peak heights show that the ME approach has a better agreement with the details of the diffraction pattern in figure 5(a).

Zoom In Zoom Out Reset image size

Compared to the dynamical method, the trajectory approaches have fewer nested loops in programming so that they can save computing time and memory when the number of Fourier expansion items is large for more accurate calculation. Besides, in equation (20) of the ME approach, there are only additions and multiplications while there are exponentiations and divisions in the other two approaches. Thus, the ME approach enables better computation efficiency. A computation time test of the three methods under the same condition is depicted in figure 9. When the number of reciprocal vector is small, the running time of the three methods is quite close. But as the size of matrices increases, the trajectory methods show significant advantages in computation efficiency. When the number of reciprocal vector grows to 21³, which is enough for accurate EBSD simulation, the Bohmian method and the ME approach take only a quarter and one tenth of the calculation time compared with the usual EBSD dynamical simulation, respectively.

Zoom In Zoom Out Reset image size