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

The Bohmian trajectory method is employed to study electron diffraction in crystalline materials. It provides a fresh understanding of the process of electron diffraction, including traveling channels of electrons and formation of diffraction patterns. By combining it with the Bloch wave method, the electron trajectories can be calculated more efficiently than the traditional wave-packet propagation algorithm. Meanwhile, we propose a momentum expectation approach which is a good approximation method with even higher computational efficiency. Both methods result in intuitive and accurate electron trajectories for the simulation of the electron backscatter diffraction (EBSD) pattern. Excellent agreement has been obtained between the simulated trajectory distributions and the experimental EBSD pattern from Mo (001) at 20 kV, where the Kikuchi patterns and higher order Laue zone rings are characterized.


Introduction
The Bohmian mechanics is an alternative interpretation of quantum mechanics, which is different from the Copenhagen interpretation and the many-worlds interpretation. The standard quantum formalism is mostly related to Copenhagen interpretation [1] for analyzing quantum phenomena from the perspective of probability and statistics. On the other hand, the Bohmian mechanics proposes a deterministic way to interpret quantum phenomena with definite particles trajectories guided by the wave function, which allows one to understand the quantum world with the familiar classical trajectory picture.
The Bohmian mechanics, first known as the pilot wave theory, was proposed by de Broglie [2,3] even before the Copenhagen interpretation and was further developed by Bohm [4,5]. It gives exactly the same predictions for all quantum phenomena as any other interpretations and is computationally friendly as has been demonstrated [6,7]. Because of its intuition, accuracy and efficiency, the Bohmian mechanics is recognized as a practical computational tool to study quantum phenomena, such as, the two-split quantum interference [8], spin system [9], identical particles [10], atom scattering [11][12][13], entangled photons [14], etc. Bohmian trajectories have already been experimentally observed [15][16][17] with weak measurement technology, and the experimental results are consistent with the calculation [18]. Thus, in recent years the Bohmian mechanics has attracted more research interest.
On the other hand, electron diffraction techniques [19,20], e.g. electron backscatter diffraction (EBSD) [21] and convergent beam electron diffraction, are significant tools for analyzing crystalline materials with electron microscopy [22]. These techniques can provide important crystallographic and phase information, such as grain size and boundary characterization, global and local texture, phase identification and separation, which is helpful to understand the properties of the materials. Though the experimental techniques for electron diffraction have been highly advanced, the formation process of the diffraction patterns is complicated and Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. deserves to be studied furthermore. It can be expected that a dynamic model of the physical process of electron diffraction will allow us to extract additional information from the experimental results. Thus, it is necessary to develop numerical simulations that can reproduce the formation process of the characteristic diffraction patterns in the form of the gorgeous Kikuchi bands. The Kikuchi pattern [23,24] is considered as being created when the electrons are scattered by different lattice planes. It exists as series of bands geometrically following the Bragg's law [25]. However, a simple geometric model could not quantitatively explain the exact intensity distribution. In order to obtain the intensity, several dynamical models of electron diffraction have been successfully developed; for example, Winkelmann et al [26] have developed a dynamical simulation method with the Bloch wave theory, and another simulation [27] was based on the multi-slice method [28].
In our work, we attempt to reproduce the electron diffraction process with quantum trajectories in a crystal based on the Bohmian mechanics. To obtain the Bohmian trajectories, time-dependent wave-packet propagation algorithms [29] are mostly applied with the use of the split operator method [30], the multi-slice method [31,32], the quantum Hamilton-Jacobi equation method [33], etc. In these methods, the timedependent Schrödinger equation is solved by calculating the evolution of the wave function over time step by step, which consumes much CPU time and computing memory. Considering the periodicity of crystalline materials, we find an efficient way to calculate the Bohmian trajectories in electron diffraction process using the Bloch wave theory. In this paper, firstly, the parallel electron beam diffraction is simulated. The wave function, the velocity field and the Bohmian trajectories will be displayed to show how a plane wave is diffracted in a crystal. Secondly, we will present an application of this approach to the simulation of the Bohmian trajectories in EBSD, where the landing positions of electrons on a screen are treated as a representation of the Bohmian trajectories. In this simulation, the EBSD pattern will be intuitively reproduced from statistical accumulation of the points on the screen. In addition, we will propose another new method called the momentum expectation (ME) approach to calculate the Bohmian trajectories. The ME approach provides a good approximation with higher efficiency to simulate the final pattern of the trajectory distribution in EBSD.

Bloch wave theory
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, where A j ( ) and k j ( ) are respectively the amplitude and the wave vector of the Bloch wave b k r , . j ( ) ( ) In the timeindependent Schrödinger equation, where U r ( ) is a periodic potential and can be expanded as a Fourier series, by substituting equations (1) and (3) into equation (2), the renowned dynamical equation of Bethe [34] is obtained as is the electron wave vector inside the crystal and k 0 is the electron wave vector in the vacuum. Then the vector k j ( ) is decomposed as, where 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 = -and C j ( ) is a column vector containing the Fourier coefficients C . j g ( ) Equation (6) can be solved by a standard method [36]. However, for high-energy electrons, it is a good approximation [36,37] 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 j g ( ) and the eigenvalues j l ( ) all become complex. To calculate the inelastic scattering processes, thermal diffuse scattering [39,40] is evaluated using the Debye-Waller factor [41][42][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.

Bohmian mechanics
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 t r, y ( )is complex so that it can be expressed in the polar form, , e x pi , , where R t r, ( )is the amplitude and S t r, ( ) is the phase. Then the time-dependent Schrödinger equation, is separated into the real part and the imaginary part [4] as The real part, equation (13), is analogous to the classical Hamilton-Jacobi equation, where S m 2 2  ( ) and U respectively stand for kinetic energy and classical potential, and Q t R mR r, 2 is interpreted as quantum potential, which is nonlocal. The particle velocity can thus be defined as, (14) for the imaginary part can be rewritten as the continuity equation, 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, 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.

ME approach
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, Based on the wave function in the form of equation (10) we have, then the ME is accordingly obtained from the integral, The averaged velocity is defined as m v p , á ñ = á ñ which describes the average electron diffraction process in the lattice. With these equations the trajectory distribution in EBSD can be evaluated.

Results and discussion
To simulate the electron diffraction in a periodic potential, we take a thin molybdenum crystal as an example without loss of generality. Molybdenum is a body-centered cubic crystal with a lattice constant of 0.31468 nm. The Debye-Waller factor is 0.23 Å 2 , which is calculated using a fourth-degree polynomial [48]. Calculation is performed for incident electrons at energy of 20 keV, and IMFP is 25.1 nm obtained from the TPP-2 M formula [45].

Diffraction of normal incident parallel electron beams
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.
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, v ,  is much larger than the perpendicular component, v , because of the high energy of the incident electrons.
To understand how the velocity field evolves over time, we now neglect v  and focus on the distribution of vertical velocity, v , 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 v . 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.
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,    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 0 , 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.
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.

Quantum Trajectories in EBSD
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.
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 2 , 10 3 , 10 4 and 10 5  trajectories), the diffraction pattern is similarly formed from a few points on the screen to a full diffraction pattern (figures 6(a)-(e)).
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 5 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).
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 5 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 7 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  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).
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 3 , 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. Figure 8. Image cross-correlation between EBSD pattern simulated by the dynamical method ( figure 5(a)) and the trajectory density distribution pattern calculated with the ME approach (figure 7(c)) (the left surface map); image cross-correlation between EBSD pattern simulated by the dynamical method ( figure 5(a)) and the trajectory density distribution pattern calculated with the Bohmian mechanics (BM) (figure 7(d)) (the right surface map).

Conclusions
The Bohmian trajectory theory is an alternative description of the quantum process from a different perspective with the Copenhagen interpretation but is more intuitive. In this paper, we have developed a Bloch wave method to study electron diffraction in a periodic potential by calculating Bohmian trajectories. In addition, we have proposed an ME approach which is more efficient to obtain the distribution of the Bohmian trajectories in EBSD simulation. These two methods consume much less computing time and memory but give the same accurate results compared with the dynamical method. From the first simulation of the normal incident parallel electron beams, the probability density is significantly related to the velocity field, i.e. a relationship exists between the evolution of the wave function and trajectories. The velocity field and the Bohmian trajectories illustrate how electrons are dynamically diffracted in crystal. It is fascinating that the electronic Talbot effect is weakened by inelastic processes and electrons will finally converge while traveling along the channels. In the second simulation, we applied the two methods to the simulation of Bohmian trajectories in EBSD, which is more complex. The formation process of EBSD pattern is demonstrated to start from a few points on a screen to a final full image with the characteristic Kikuchi bands and HOLZ rings, which is similar to a single-electron experiment. Bohmian mechanics and the ME approach thus allow us to have a deeper understanding as well as to acquire quantitative information of electron diffraction from a trajectory perspective rather than from probability density. Technically, electron, photon and proton diffraction techniques can all be studied with these methods. Besides the application in the simulation of EBSD, it is expected to contribute to the study of crystallographic defects by the abnormal propagation of the trajectories.