## 3 Band Structure Governed Wave Packet Dynamics  
In the traditional formulation of the WPD, the environment of the moving electron quasiparticle (i.e., the crystal and the averaged effect of the other electrons) is incorporated into the potential energy operator, often by the help of a (pseudo)potential, thus the Hamiltonian is \(\hat{H}=\hat{K}_{free}+\hat{V}\), where \(\hat{K}_{free}\) is the kinetic energy operator for a free electron. In this section, we present an alternative formulation, where the effect of the crystal potential and the many-body effects are entirely incorporated into the kinetic energy operator, i.e., the Hamiltonian is \(\hat{H}=\hat{K}+\hat{V}_{free}\). In our alternative formulation all crystal structure and electronic structure effects are taken care of in band structure calculations that yield a dispersion relation \(E_{n}\left(\overset{\rightarrow}{k}_{Bloch}\right)\), where \(\overset{\rightarrow}{k}_{Bloch}\) is the Bloch wave vector and \(n\) is the band index. Thus, we replace the kinetic energy operator in free space, \(\hat{K}_{free}=\dfrac{\left|\overset{\rightarrow}{k}\right|^{2}}{2}=E_{free} \left(\overset{\rightarrow}{k}\right)\) in momentum representation, where \(E_{free}\left(\overset{\rightarrow}{k}\right)\) is the free space dispersion relation, with the kinetic energy operator that describes the dispersion of a many-electron system in a given crystalline material. Thus, the momentum representation of the kinetic energy operator becomes \(\hat{K}=E\left(\overset{\rightarrow}{k}_{Bloch}\right)\). By using this modified kinetic operator \(\hat{K}=E\left(\overset{\rightarrow}{k}_{Bloch}\right)\), a perfect crystal Hamiltonian (without any defect) can be written as \(\hat{H}=\hat{K}+\hat{V}_{free}\), where \(\hat{V}_{free}=0\) by definition. We have already presented such calculations for perfect 2D crystals, graphene and MoS\({}_{2}\) single sheets in Reference [20]. Those calculations successfully reproduced the trigonal warping effect and the anisotropic WP spreading characteristic of these 2D materials and showed different symmetries of the WP spreading, depending on the band structure and the spectral distribution of the initial WP. As we emphasized in Reference [20], similar calculations can be easily performed for any crystalline material, where the dispersion relation is known. In this section, we combine the two methods to study defects in 2D crystalline materials. We describe the (infinite, periodic) crystal by the kinetic energy operator (utilizing the dispersion relation of the material) and the structural defect by a local potential. Thus,  
Figure 3: Simulated and measured FFT-LDOS maps. (**a**) Two-dimensional Fourier transform (2D FFT) of the probability density of the initial wave packet shown on Figure 2d. (**b**) 2D-FFT of the probability density at t = 7.2 fs (Figure 2f). (**c**) 2D-FFT of a topographic STM image. our Hamiltonian will become \(\hat{H}=\hat{K}_{crystal}+\hat{V}_{defect}\), where the \(\hat{K}_{crystal}=E\bigg{(}\overset{\rightarrow}{k}_{Blch}\bigg{)}\) kinetic energy operator describes the electronic structure of the crystal lattice and the \(\hat{V}_{defect}=V\bigg{(}\overset{\rightarrow}{r}\bigg{)}\) potential energy operator describes the local defect. The main novelty of this combined method is the straightforward treatment of the different types of defects in several 2D materials. Figure 4a-c show the time evolution of a WP for an unperturbed graphene lattice and Figure 4d-f show the time evolution in the presence of a local defect calculated by numerically solving the time-dependent Schrodinger equation. The defect was modeled by a Gaussian potential, its height is 2 eV and its HWHM is 0.8 nm. In contrast to the previous Section, now we used a localized initial state. The spectral distribution of the initial WP is a sum of six Gaussians, placed in the six \(K\) and \(K^{\prime}\) points (Figure 4g) resulting in a spatially well-localized initial state.