Dynamically Characterizing the Structures of Dirac Points via Wave Packets

Dan-Dan Liang(梁丹丹)$^{1,2}$, Xin Shen(沈鑫)$^{3\hyperlink{s}{*}}$, and Zhi Li(李志)$^{1,2\hyperlink{s}{*}}$

doi:10.1088/0256-307X/40/11/110302

⇦Chinese Physics Letters, 2023, Vol. 40, No. 11, Article code 110302⇨ Dynamically Characterizing the Structures of Dirac Points via Wave Packets Dan-Dan Liang (梁丹丹)1,2, Xin Shen (沈鑫)3*, and Zhi Li (李志)1,2* Affiliations 1Key Laboratory of Atomic and Subatomic Structure and Quantum Control (Ministry of Education), Guangdong Basic Research Center of Excellence for Structure and Fundamental Interactions of Matter, School of Physics, South China Normal University, Guangzhou 510006, China 2Guangdong Provincial Key Laboratory of Quantum Engineering and Quantum Materials, Guangdong-Hong Kong Joint Laboratory of Quantum Matter, Frontier Research Institute for Physics, South China Normal University, Guangzhou 510006, China 3College of Sciences, China Jiliang University, Hangzhou 310018, China Received 6 July 2023; accepted manuscript online 25 September 2023; published online 18 October 2023 ^*Corresponding authors. Email: shenx@cjlu.edu.cn; lizphys@m.scnu.edu.cn Citation Text: Liang D D, Shen X, and Li Z 2023 Chin. Phys. Lett. 40 110302 Abstract Topological non-trivial band structures are the core problem in the field of topological materials. We investigate the topological band structure in a system with controllable Dirac points from the perspective of wave packet dynamics. By adding a third-nearest-neighboring coupling to the graphene model, additional pairs of Dirac points emerge. The emergence and annihilation of Dirac points result in hybrid and parabolic points, and we show that these band structures can be revealed by the dynamical behaviors of wave packets. In particular, for the gapped hybrid point, the motion of the wave packet shows a one-dimensional Zitterbewegung motion. Furthermore, we also show that the winding number associated with the Dirac point and parabolic point can be determined via the center of mass and spin texture of wave packets, respectively. The results of this work could motivate new experimental methods to characterize a system's topological signatures through wave packet dynamics, which may also find applications in systems of other exotic topological materials.

DOI:10.1088/0256-307X/40/11/110302 © 2023 Chinese Physics Society Article Text Ever since the birth of quantum mechanics, the solid band theory has provided a fundamental explanation for metals, insulators and semiconductors.[1] Nevertheless, both theoretical and experimental progresses in condensed matter physics have greatly extended the classification of materials and deepened the concept of phase transition in recent years.[2,3] Among the novel physics are topological insulators and superconductors[4,5] or semimetals,[6-8] which have inspired many investigations in quantum transport[9] and synthetic topological quantum matter,[10-14] and hence potential applications in quantum computation.[15] Among topological materials, the band structure of a Dirac point is ubiquitous in indicating the topological phase transition or band inversion.[16,17] A band's topology can be characterized by an invariant.[18] Much effort has been paid to exploring the topological states via measurements of topological invariants or edge states, especially in the experimental simulation of topological matters in various platforms such as optical crystals,[12] trapped ions,[14] and ultracold atoms.[13] In this Letter, we report the investigation of band structure via wave packet dynamics that uncovers the topological property encoded in the band. The corresponding Dirac points can be simulated in an optical lattice in the form of low-energy excitation.[19] The number of Dirac points is restricted to even numbers with opposite winding numbers or chirality[20] due to the Fermi doubling theorem.[21] By tuning the laser intensity, the motion of Dirac points can be induced with the emergence or annihilation of Dirac points. We focus on a prototypical graphene model hosting Dirac quasiparticles with third-nearest-neighboring (N3) coupling considered.[22,23] The N3 coupling shifts the position of the Dirac points and modifies the band structure into a topologically trivial hybrid point or a topologically non-trivial parabolic point with a higher winding number.[24] A wave packet can be viewed as a state with a definite (quasi)momentum, therefore its expansion reveals local band structure such as the group velocity and band gap. The merging of Dirac points with opposite winding numbers gives rise to a gapped hybrid point, leading to a directed Zitterbewegung (ZB) motion.[25] The merging process also happens for Dirac points with the same winding number, forming a gapless parabolic point of a higher winding number. From the dynamics of a single wave packet we show that the winding number can be inferred from the motion of the wave packet's center of mass and spin texture, providing a dynamical framework to characterize the topological band structures.[26-29] Model and Band Structures. Within the tight-binding approximation, the Hamiltonian of the graphene can be written as $H=\sum_{mn}c^†_m J_{mn} c_n$, where $c$ and $c^†$ are the creation and annihilation operators of each lattice site, and $J_{mn}$ characterize the hopping between lattices sites. Considering only the nearest neighboring hopping, gapless Dirac points emerge in the graphene merely at the ${\boldsymbol K}$ and ${\boldsymbol K'}$ valleys. By adding N3 coupling to the tight-binding model of graphene as shown in Fig. 1, the preserved chiral symmetry allows for the emergence of additional pairs of Dirac points. The corresponding Hamiltonian reads \begin{align} H=J\sum_{i,j=NN}c_{A,i}^† c_{B,j}+J_3\sum_{i,j=N3}c_{A,i}^† c_{B,j}+{\rm h.c.}, \tag {1} \end{align} where $\rm{NN}$ stands for the hopping between the nearest neighboring sites and $\rm{N3}$ characterizes the hopping between the N3 sites with hopping parameters $J$ and $J_3$, respectively. $A$ and $B$ denote two sublattices in the monolayer graphene. By Fourier transformation, the corresponding Bloch Hamiltonian reads \begin{align} \mathcal H_{\boldsymbol k}= \begin{bmatrix} 0 & f_{\boldsymbol k}\\ f_{\boldsymbol k}^* & 0 \end{bmatrix}, \tag {2} \end{align} where \begin{align} f_{\boldsymbol{k}}=\,&J\big(1+e^{i\boldsymbol{k}\cdot \boldsymbol{a}_1}+e^{i\boldsymbol{k}\cdot \boldsymbol{a}_2}\big)\notag\\ &+J_3\big[e^{i\boldsymbol{k}\cdot (\boldsymbol{a}_1+ \boldsymbol{a}_2)}+e^{i\boldsymbol{k}\cdot (\boldsymbol{a}_1- \boldsymbol{a}_2)}+e^{i\boldsymbol{k}\cdot (\boldsymbol{a}_2-\boldsymbol{a}_1)}\big], \tag {3} \end{align} $\boldsymbol{a}_1=(\frac{\sqrt{3}a}{2}, \frac{3a}{2})$ and $\boldsymbol{a}_2=(-\frac{\sqrt{3}a}{2}, \frac{3a}{2})$ are the elementary vectors of the lattice with lattice constant $a$; ${\boldsymbol k}$ is the quasimomentum. Due to the $J_3$ term, the energy spectrum is drastically changed. When $J_3$ increases and reaches the critical value $J/3$, a new pair of Dirac points emerge from each of the three inequivalent time-reversal invariant $\boldsymbol{M}$ points in reciprocal space. When $J_3$ is further increased either of the newly emerged pair of Dirac points moves the ${\boldsymbol M}$ point to the ${\boldsymbol K}$ or ${\boldsymbol K'}$ valley and then four Dirac points merge into a single degeneracy. The critical value is $J_3=J/2$, in which case the dispersion becomes quadratic. Near $J_3\approx{J/2}$, the Hamiltonian in the vicinity of the $\boldsymbol{K'}$ point can be derived as \begin{align} \mathcal H_K= \begin{bmatrix} 0 & -\frac{\boldsymbol{q}^2}{2m^*}+c\boldsymbol{q}^†+\varDelta\\ -\frac{\boldsymbol{q}^{†2}}{2m^*}+c\boldsymbol{q}+\varDelta^* & 0 \end{bmatrix}, \tag {4} \end{align} where ${\boldsymbol q} \equiv q_x+iq_y$ and ${\boldsymbol k}={\boldsymbol K'}+{\boldsymbol q}$. Here, ${\boldsymbol q}$ denotes the displacement from ${\boldsymbol K'}$ point. The parameters are $m^*=-4/9J$, $c=3J/2-3J_3$ and $\varDelta=0$. When $J_3=J/2$, i.e., $c=0$, the dispersion shows quadratic behavior in both momenta $q_x$ and $q_y$. This differs from the usual parabolic kinetic energy and carries a vorticity.

Fig. 1. (a) The tight-binding model: schematic diagram of a honeycomb lattice with nearest and third-nearest neighboring hopping. (b) Brillouin zone and the high symmetry points.

Fig. 2. Band structures (left panels) of Hamiltonian (4) and the corresponding density profiles $|\psi({\boldsymbol r},t)|^2$ (right panels) of the wave packets, with $\varDelta=-2\varDelta_{\rm c}$, $-\varDelta_{\rm c}$, $-0.5\varDelta_{\rm c}$, $0$, in (a)–(d) and (e)–(h). The red point marks the center of the wave packet. The units for momentum, position and time are $1/a$, $a$, and $a/c$, respectively. The width of the wave packets is $d=5a$. Throughout the simulation $m=1/ac$.

In the following, the parameter $\varDelta$ is assigned with a purely imaginary value $\varDelta=ib$ ($b$ is a real number), such that all the band characteristics are included in the Hamiltonian given in Eq. (4), i.e., emergence or annihilation of Dirac points with opposite chirality, and also a quadratic dispersion carrying nontrivial vorticity.[24] The energy spectra are plotted in Fig. 2. Starting from Fig. 2(d), when $\varDelta = 0$ there are four Dirac cones. Increasing or decreasing $\varDelta$ causes two Dirac points to merge, resulting in a linear–quadratic degeneracy and a gapped band, also known as a hybrid point. This can be seen from Figs. 2(a)–2(d). Next, we will discuss the band structures of the model for varying values of $\varDelta$, along with their corresponding dynamical behaviors. The Wave Packet Dynamics. The wave function for the low-energy quasiparticles satisfies the Schrödinger's equation $i \partial_t \psi={\mathcal H}_{K}\psi$ ($\hbar=1$). Throughout the dynamical simulation, we take the initial state as a two-dimensional Gaussian wave packet.[30-32] The wave function in momentum space is \begin{align} \boldsymbol{\psi}_{\boldsymbol k}(t=0)=d\sqrt{\frac{2}{\pi}} e^{-d^2[(k_x-k_1)^2+(k_y-k_2)^2]}{\phi}, \tag {5} \end{align} where $d$ is the width of the wave packet and $k_1(k_2)$ is the initial momentum in the $x(y)$ direction; $\phi=(c_1,c_2)^{\scriptscriptstyle{\rm T}}$ is the spinor part with superscript T standing for the matrix transposition. When the width of a wave packet $d$ is sufficiently large, it is narrowly spread around the point $(k_1,k_2)$ in momentum space, and its time evolution is governed by the Hamiltonian near that point. In order to reveal the local Dirac band structures of Eq. (4), $(k_1,k_2)$ take values of the Dirac points or hybrid points marked in Fig. 2. For a specific demonstration, a monocomponent spinor $\phi=(1,0)^{\scriptscriptstyle{\rm T}}$ is considered. To understand the motion behavior, we proceed to calculate the expectation value of the position operator. By expressing the Hamiltonian (4) as $\mathcal H_K=d_x\sigma_x+d_y\sigma_y$, the position of the center of mass (PCM) is obtained as \begin{align} \bar{x}(t)=\,&\langle \psi(t)|\hat x|\psi(t)\rangle\notag\\ =\,&\frac{2d^2}{\pi}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\frac{\frac{q_y}{m^*}d_x-(c-\frac{q_x}{m^*})d_y}{2E^2}[\cos(2\,Et)-1]\notag\\ &\times e^{-2d^2[(k_x-k_1)^2+(k_y-k_2)^2]}dk_xdk_y,\notag\\ \bar{y}(t)=\,&\langle \psi(t)|\hat y|\psi(t)\rangle\notag\\ =\,&\frac{2d^2}{\pi}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\frac{(c+\frac{q_x}{m^*})d_x-\frac{q_y}{m^*}d_y}{2E^2}[\cos(2\,Et)-1]\notag\\ &\times e^{-2d^2[(k_x-k_1)^2+(k_y-k_2)^2]}dk_xdk_y, \tag {6} \end{align} where $E=\sqrt{d^2_x+d^2_y}$. The PCM can also be obtained by directly solving Schrödinger's equation, and the wave packet's density profile versus time is shown in Fig. 2. Both the analytical and numerical results for the PCM are plotted in Fig. 3. From these diagrams, the local band structures can be inferred.

Fig. 3. Analytically (lines) and numerically (symbols) calculated PCM $\bar x$ (a) and $\bar y$ (b) with different values of $\varDelta$. The width of the initial wave packet $d=5a$.

As shown in Fig. 2(e), the wave packet exhibits oscillatory motion, corresponding to the gapped band structure after the annihilation of Dirac points. Unlike the two-dimensional massive Dirac quasiparticles, the wave packet here oscillates one-dimensionally due to the linear–quadratic character of the dispersion. When the gap closes, the oscillation ceases and the wave packet splits into two sub-packets moving in opposite directions [see Figs. 2(b) and 2(f)]. In this scenario, the one-dimensional motion results from a finite group velocity in the linear direction but zero in the quadratic direction. In Fig. 2(c), the wave packet shows diffusion due to finite velocities in all directions, corresponding to the band structure of the newly emerged Dirac point. The density distribution is anisotropic because the group velocity in the previous quadratic direction is relatively smaller than the velocity in the linear direction. When $\varDelta=0$, the four Dirac points are located in momentum space with $C_3$ symmetry, and the initial isotropic wave packet now features $C_3$ symmetry as well [see Fig. 2(d)]. The smaller the width, the broader the wave packet in momentum space, resulting in the inclusion of more momenta in the dynamics. The dynamics of the local band structure would then be smeared out for a narrow wave packet in position space. In Fig. 3 with $\varDelta = -2\varDelta_{\rm c}$, the oscillation shows damping behavior due to the finite width of the wave packets. The dynamics can be modeled as a group of oscillators, with each frequency corresponding to a ${\boldsymbol k}$-dependent energy gap. To manifest the dynamics of the local band structure, a proper width can be inferred from the ZB dynamics with $d\gg c/\varDelta\sim a$.[33,34] In the optical lattice, the trapping frequency is typically of the order $10^2$ Hz, and the lattice constant $\sim$ $\lambda/2$, with $\lambda$ being the laser wavelength. For the trapped $^{87}$Rb Bose condensate,[35] for example, the ratio $d/a\sim 17$ is much larger than that in the numerical simulations shown in Figs. 2 and 3. PCM and Winding Number. By tuning the parameter $\varDelta$, the merging process involves a phase transition from a semimetal to an insulator, where two massless Dirac points of opposite chirality combine into a topologically trivial gapped point. In its most general form, a gapless Dirac Hamiltonian has the form \begin{align} \mathcal H_D({\boldsymbol q})={\boldsymbol q}\cdot({\boldsymbol v}_1\sigma_x+{\boldsymbol v}_2\sigma_y). \tag {7} \end{align} The Berry connection is defined as $A_i=\langle u_{\boldsymbol p} | \partial_i |u_{\boldsymbol p}\rangle$, where $|u\rangle$ is Hamiltonian's eigenstate. For a Dirac point, the winding number is $w=\oint {\boldsymbol A}\cdot d{\boldsymbol l}$, in which the integration loop circles around the degenerate point. For the Hamiltonian (7), the winding is derived as $w=\lambda{\rm sgn}({\boldsymbol v}_1\wedge{\boldsymbol v}_2)$, where $\lambda=\pm 1$ is the band index and the wedge product here is defined as ${\boldsymbol v}_1 \wedge {\boldsymbol v}_2\equiv ({\boldsymbol v}_1\times {\boldsymbol v}_2)_z$. Therefore, we have the relationship between the winding number and the velocities ($\lambda=1$) as follows: \begin{align} w={\rm sgn}({\boldsymbol v}_1\wedge{\boldsymbol v}_2)= \left \{ \begin{array}{ll} \!+1, ~~ {\rm for} ~ v_{1x}v_{2y} > v_{1y}v_{2x},\\ \!-1, ~~ {\rm for} ~ v_{1x}v_{2y} < v_{1y}v_{2x}. \end{array} \right. \tag {8} \end{align} Similarly, we consider the evolution of the wave packet in the form of Eq. (5) with $\phi=(c_1,c_2)^{\scriptscriptstyle{\rm T}}$ and $k_1=k_2=0$. In the large width limit, the spatial dynamics is totally determined by the Hamiltonian ${\mathcal H}_D({\boldsymbol q}=0)$. The PCM accordingly can be obtained as \begin{align} \bar{x}(t)&=[(c_1c_2^*+c_1^*c_2)v_{1x}+i(c_1c_2^*-c_1^*c_2)v_{2x}]t,\notag\\ \bar{y}(t)&=[(c_1c_2^*+c_1^*c_2)v_{1y}+i(c_1c_2^*-c_1^*c_2)v_{2y}]t. \tag {9} \end{align} One can then find two initial spinor states, for example, $|\phi_1 \rangle=\frac{1}{\sqrt{2}}(1,1)^{\scriptscriptstyle{\rm T}}$ and $|\phi_2 \rangle=\frac{1}{\sqrt{2}}(1,i)^{\scriptscriptstyle{\rm T}}$, such that the corresponding PCMs of the wave packets take a simple form \begin{align} \tag {10} &\bar{x}_1=v_{1x}t, ~~~\bar{y}_1=v_{1y}t,\\ & \tag {11} \bar{x}_2=v_{2x}t, ~~~\bar{y}_2=v_{2y}t, \end{align} respectively. In a straightforward manner, we rewrite the winding number in terms of PCM as \begin{align} w={\rm sgn}(\bar{x}_1\bar{y}_2-\bar{x}_2\bar{y}_1)= \left \{ \begin{array}{ll} \!+1, ~~~ {\rm for}~ \bar{x}_1\bar{y}_2>\bar{x}_2\bar{y}_1,\\ \!-1, ~~~ {\rm for}~ \bar{x}_1\bar{y}_2 < \bar{x}_2\bar{y}_1. \end{array} \right. \tag {12} \end{align} Numerically, we apply the above method to demonstrate the winding number of the four Dirac points when $\varDelta=0$, in which case four different Dirac points are located at $(0,0)$, $(2m^*c,0)$, $(-m^*c,-\sqrt{3}m^*c)$ and $(-m^*c,\sqrt{3}m^*c)$. Two of them are plotted in Fig. 4, from which the winding number can be directly inferred.

Fig. 4. Analytically (lines) and numerically (symbols) calculated $\bar x\bar y$ as a function of time $t$. Here $(\bar x_1, \bar y_1)$ and $(\bar x_2, \bar y_2)$ are calculated with initial states $|\phi_1 \rangle=1/\sqrt{2}(1,1)^{\scriptscriptstyle{\rm T}}$ and $|\phi_2 \rangle=1/\sqrt{2}(1,i)^{\scriptscriptstyle{\rm T}}$, respectively. For (a) and (b), the momentum center of the wave packet is located at the Dirac points $(0,0)$, $(2m^*c,0)$, respectively. The width of the initial wave packet $d=15a$.

In the case that ${\boldsymbol v}_1=(v_{1x},0)$ and ${\boldsymbol v}_2=(0,v_{2y})$, which is commonly seen in topological insulators, via Eq. (9) the PCM is now \begin{align} \bar x(t)=&2R\cos\theta \,v_{1x}t,\notag\\ \bar y(t)=&2R\sin\theta \,v_{2y}t, \tag {13} \end{align} where $c_1^*c_2 \equiv Re^{i\theta}$. The winding number around the Dirac point is simply \begin{align} w={\rm sgn}(v_{1x}v_{2y})= \left\{ \begin{array}{ll} \!+1, ~~~ {\rm for}~ v_{1x}v_{2y}>0,\\ \!-1, ~~~ {\rm for}~ v_{1x}v_{2y} < 0. \end{array} \right. \tag {14} \end{align} Thus, we can determine the sign of $(v_{1x},v_{2y})$ and obtain the winding number by analyzing the motion direction of the PCM with respect to the initial state parameter $\theta$. Spin Texture and Winding Number. Within topological band theory, the winding number can also be viewed as a mapping between manifolds. Rewriting the Hamiltonian (7) in the form of ${\mathcal H}_D=h{\boldsymbol n}({\boldsymbol k})\cdot {\boldsymbol \sigma}$, the associated winding number is[20] \begin{align} w=\frac{1}{\pi}\oint_\mathcal C n_x d n_y, \tag {15} \end{align} where $\oint_\mathcal C$ denotes the integration along the loop circling around the Dirac point. The Bloch vector lies along the equator and the winding number $w$ characterizes the mapping from the momentum loop to the wave function space ${\boldsymbol n}$. The Hamiltonian can be seen as a magnetic field acting on a spin while its direction ${\boldsymbol n}$ varies in momentum space. The number of times that the vector ${\boldsymbol n}$ winds along the equator gives the winding number $w$. Here the initial state wave function is a Gaussian wave packet with a finite width and a spin-up spinor throughout momentum space. The evolution of the wave packet is the collection of spin precessions around ${\boldsymbol n}$ at each momentum. The initial spinor state is perpendicular to the ${\boldsymbol n}$, hence the spin will always be perpendicular to it during the evolution. By examining the spin texture of the wave packet on isoenergy surfaces, one can derive the winding of ${\boldsymbol n}$ from the spin texture information. For example, considering the case of $\varDelta=0$ and $c=0$ in Eq. (4), the corresponding Hamiltonian reads \begin{align} \mathcal H_V= \begin{bmatrix} 0 &-\frac{q^2}{2m^*}e^{i2\varphi}\\ -\frac{q^2}{2m^*}e^{-i2\varphi} & 0 \end{bmatrix}, \tag {16} \end{align} where $q^2=q_x^2+q_y^2$ and $\tan \varphi=q_y/q_x$. The spin texture, ${\boldsymbol s}\equiv\langle \psi_{\boldsymbol q} | {\boldsymbol \sigma}(t)| \psi_{\boldsymbol q} \rangle$, projected in the $xy$-plane can be obtained as \begin{align} s_x({\boldsymbol q},t)&=\sin(2\varphi)\sin(2Et),\notag\\ s_y({\boldsymbol q},t)&=\cos(2\varphi)\sin(2Et). \tag {17} \end{align} In Fig. 5, we plot the spin texture of both the Dirac point and the parabolic point. As only the direction of the spin texture determines the winding number, the overall time-dependent factor has been ignored.

Fig. 5. The spin texture for (a) the Dirac point and (b) the parabolic point with winding numbers $w=1$ and $w=2$, respectively. The width of wave packet is $d=5a$.

Experimental Implementation. Finally, we discuss the experimental implementation of initial state preparation and detection. In the cold atom system, the Bose gas confined in a harmonic trap can be approximately modeled by a Gaussian function, which naturally serves as a wave packet for simulation purposes here. By applying a magnetic field gradient[19] or moving the optical lattice,[36,37] the atomic cloud can be accelerated such that the quasimomentum increases linearly up to the expected value. After releasing the trap and switching on the optical lattice potential, the atomic ensemble would start expanding in the two-dimensional honeycomb lattice. The time-evolving density profile can then be measured by the absorption imaging technique, hence allowing for the measurement of PCM.[38] For the measurement of the spin texture, the distribution ${\boldsymbol s}({\boldsymbol q})$ can be experimentally determined from the time-of-flight images after the optical lattice is switched off, through which the momentum density distribution can be obtained. By the Raman impulse, $s_x({\boldsymbol q})$ and $s_y({\boldsymbol q})$ can be mapped to $s_z({\boldsymbol q})$, which is directly related to the (pseudo)spin-up and (pseudo)spin-down density $s_z({\boldsymbol q})=\frac 1 2[n_{\uparrow}({\boldsymbol q})+n_{\downarrow}({\boldsymbol q})]/[n_{\uparrow}({\boldsymbol q})+n_{\downarrow}({\boldsymbol q})]$.[39-41] Preparation of the initial spinor state can be easily achieved by utilizing light–atom interaction if the atomic internal state is associated with the pseudospin degree of freedom. In the graphene model here, the spin refers to the sublattice. When the width of the wave packet is much larger than the lattice constant, the neighboring sites are equally populated and thus the initial spinor takes the form $(1,1)$.[42] For any other spinor state, it can be prepared by switching on the optical lattice and driving the atoms to a proper quasimomentum. By using the lattice Hamiltonian near the quasimomentum, the desired spinor can then be prepared for a specific holding time. Since now the wave packet engages with only the spatial degrees of freedom, the spin texture of the wave packet can then be resolved through the band tomography method.[43,44] In summary, we have investigated the wave packet dynamics of the tight-binding graphene model, which includes third-nearest-neighboring coupling. The motion of Dirac points leads to band structures with gapped or gapless hybrid points. Corresponding wave packet dynamics show directed ZB motion and $C_3$ symmetry, indicating linear–quadratic and $C_3$ symmetric dispersion, respectively. We have also built up a relationship between the winding number of the Dirac point and the motion of the wave packet, and show that the winding numbers of both the Dirac point and parabolic points can be characterized by the spin texture of the wave packet. The elucidation of wave packet dynamics has been extensively performed for the simulation of relativistic phenomena in various quantum and classical setups. Therefore, we could expect that our results will benefit the experimental simulation of topological materials. Moreover, new topological states of matter, such as the higher-order topological insulators and semimetals,[45] topological Euler insulators[11,46] and non-Hermitian topological insulators in open systems, have recently attracted increasing interest.[47] Due to the similarity or novel features of the band structure, our results also provide a dynamical perspective for the exploration of the Dirac points[48] or topological characteristics in those topological quantum states.[49,50] Acknowledgements. This work was supported by the National Key Research and Development Program of China (Grant No. 2022YFA1405300), the National Natural Science Foundation of China (Grant Nos. 12074180 and 12104430), and the Guangdong Basic and Applied Basic Research Foundation (Grant No. 2021A1515012350). References Colloquium : Zoo of quantum-topological phases of matterColloquium : Topological band theoryColloquium : Topological insulatorsTopological insulators and superconductorsDirac Semimetal in Three DimensionsDirac Semimetals in Two DimensionsTopological semimetal and Fermi-arc surface states in the electronic structure of pyrochlore iridatesQuantum transport in topological semimetals under magnetic fieldsTopological AcousticsExperimental characterization of fragile topology in an acoustic metamaterialTopological photonicsTopological quantum matter with cold atomsProgrammable quantum simulations of spin systems with trapped ionsNon-Abelian anyons and topological quantum computationTopological characterization of a one-dimensional optical lattice with a forceLink between Zitterbewegung and topological phase transitionsClassification of topological quantum matter with symmetriesCreating, moving and merging Dirac points with a Fermi gas in a tunable honeycomb latticeAn Aharonov-Bohm interferometer for determining Bloch band topologyThe Adler-Bell-Jackiw anomaly and Weyl fermions in a crystalDirac point metamorphosis from third-neighbor couplings in graphene and related materialsAn equivalence between monolayer and bilayer honeycomb latticesWinding Vector: How to Annihilate Two Dirac Points with the Same ChargeTopological bands for ultracold atomsSynthetic Spin-Orbit Coupling in Two-Level Cold AtomsProbe Knots and Hopf Insulators with Ultracold AtomsPhase-Modulated 2D Topological Physics in a One-Dimensional Ultracold SystemZitterbewegung for ultracold atoms in the merging of Dirac pointsDynamics of Weyl quasiparticles in an optical latticeSimulating the Majorana dynamics with ultracold atomic gases in a bilayer honeycomb latticeAtomic ZitterbewegungZitterbewegung (trembling motion) of electrons in semiconductors: a reviewAtomic Bose–Einstein condensate in twisted-bilayer optical latticesObservation of Dynamical Instability for a Bose-Einstein Condensate in a Moving 1D Optical LatticeAtomic Bloch-Zener Oscillations and Stückelberg Interferometry in Optical LatticesDirect observation of zitterbewegung in a Bose–Einstein condensateSeeing Topological Order in Time-of-Flight MeasurementsRealization of two-dimensional spin-orbit coupling for Bose-Einstein condensatesWave Packet Dynamics in Synthetic Non-Abelian Gauge FieldsLandau-Zener-Stückelberg interferometry in

PT

-symmetric non-Hermitian modelsExperimental reconstruction of the Berry curvature in a Floquet Bloch bandBloch state tomography using Wilson linesHigher-order band topologyFragile Topology and Wannier ObstructionsNon-Hermitian Topological Phenomena: A ReviewTunable Dirac points in a two-dimensional non-symmorphic wallpaper group latticeUncover Topology by Quantum Quench DynamicsScheme to Measure the Topological Number of a Chern Insulator from Quench Dynamics

[1]	Ashcroft N W and Mermin N D 1976 Solid State Physics (Saunders College Publishing)
[2]	Wen X G 2017 Rev. Mod. Phys. 89 041004
[3]	Bansil A, Lin H, and Das T 2016 Rev. Mod. Phys. 88 021004
[4]	Hasan M Z and Kane C L 2010 Rev. Mod. Phys. 82 3045
[5]	Qi X L and Zhang S C 2011 Rev. Mod. Phys. 83 1057
[6]	Young S M, Zaheer S, Teo J C Y, Kane C L, Mele E J, and Rappe A M 2012 Phys. Rev. Lett. 108 140405
[7]	Young S M and Kane C L 2015 Phys. Rev. Lett. 115 126803
[8]	Wan X G, Turner A M, Vishwanath A, and Savrasov S Y 2011 Phys. Rev. B 83 205101
[9]	Lu H Z and Shen S Q 2017 Front. Phys. 12 127201
[10]	Yang Z J, Gao F, Shi X H, Lin X, Gao Z, Chong Y D, and Zhang B L 2015 Phys. Rev. Lett. 114 114301
[11]	Peri V, Song Z D, Serra-Garcia M, Engeler P, Queiroz R, Huang X, Deng W, Liu Z, Bernevig B A, and Huber S D 2020 Science 367 797
[12]	Ozawa T, Price H M, Amo A, Goldman N, Hafezi M, Lu L, Rechtsman M C, Schuster D, Simon J, Zilberberg O, and Carusotto I 2019 Rev. Mod. Phys. 91 015006
[13]	Zhang D W, Zhu Y Q, Zhao Y X, Yan H, and Zhu S L 2018 Adv. Phys. 67 253
[14]	Monroe C, Campbell W C, Duan L M, Gong Z X, Gorshkov A V, Hess P W, Islam R, Kim K, Linke N M, Pagano G, Richerme P, Senko C, and Yao N Y 2021 Rev. Mod. Phys. 93 025001
[15]	Nayak C, Simon S H, Stern A, Freedman M, and Sarma S D 2008 Rev. Mod. Phys. 80 1083
[16]	Shen X and Li Z 2018 Phys. Rev. A 97 013608
[17]	Shen X, Zhu Y Q, and Li Z 2022 Phys. Rev. B 106 L180301
[18]	Chiu C K, Teo J C Y, Schnyder A P, and Ryu S 2016 Rev. Mod. Phys. 88 035005
[19]	Tarruell L, Greif D, Uehlinger T, Jotzu G, and Esslinger T 2012 Nature 483 302
[20]	Duca L, Li T, Reitter M, Bloch I, Schleier-Smith M, and Schneider U 2015 Science 347 288
[21]	Nielsen H and Ninomiya M 1983 Phys. Lett. B 130 389
[22]	Bena C and Simon L 2011 Phys. Rev. B 83 115404
[23]	Montambaux G 2012 Eur. Phys. J. B 85 375
[24]	Montambaux G, Lim L K, Fuchs J N, and Piéchon F 2018 Phys. Rev. Lett. 121 256402
[25]	Schrödinger E 1930 Sitzungsber. Preuss. Akad. Wiss. Phys. Math. Kl. 24 418
[26]	Cooper N R, Dalibard J, and Spielman I B 2019 Rev. Mod. Phys. 91 015005
[27]	Zhang Q, Gong J B, and Oh C H 2013 Chin. Phys. Lett. 30 080301
[28]	Deng D L, Wang S T, Sun K, and Duan L M 2018 Chin. Phys. Lett. 35 013701
[29]	Guo G F, Bao X X, Tan L, and Gu H Q 2021 Chin. Phys. Lett. 38 040302
[30]	Li Z, Cao H, and Fu L B 2015 Phys. Rev. A 91 023623
[31]	Li Z, Wang H Q, Zhang D W, Zhu S L, and Xing D Y 2016 Phys. Rev. A 94 043617
[32]	Shen X, Zhang D W, Yan H, Li Z, and Zhu S L 2020 Phys. Rev. Res. 2 013037
[33]	Merkl M, Zimmer F E, Juzeliūnas G, and öhberg P 2008 Europhys. Lett. 83 54002
[34]	Zawadzki W and Rusin T M 2011 J. Phys.: Condens. Matter 23 143201
[35]	Meng Z M, Wang L W, Han W, Liu F D, Wen K, Gao C, Wang P J, Chin C, and Zhang J 2023 Nature 615 231
[36]	Fallani L, Sarlo L D, Lye J E, Modugno M, Saers R, Fort C, and Inguscio M 2004 Phys. Rev. Lett. 93 140406
[37]	Kling S, Salger T, Grossert C, and Weitz M 2010 Phys. Rev. Lett. 105 215301
[38]	LeBlanc L J, Beeler M C, Jiménez-Garcĺa K, Perry A R, Sugawa S, Williams R A, and Spielman I B 2013 New J. Phys. 15 073011
[39]	Alba E, Fernandez-Gonzalvo X, Mur-Petit J, Pachos J K, and Garcia-Ripoll J J 2011 Phys. Rev. Lett. 107 235301
[40]	Wu Z, Zhang L, Sun W, Xu X T, Wang B Z, Ji S C, Deng Y, Chen S, Liu X J, and Pan J W 2016 Science 354 83
[41]	Hasan M, Madasu C S, Rathod K D, Kwong C C, Miniatura C, Chevy F, and Wilkowski D 2022 Phys. Rev. Lett. 129 130402
[42]	Shen X, Wang F, Li Z, and Wu Z 2019 Phys. Rev. A 100 062514
[43]	Fläschner N, Rem B S, Tarnowski M, Vogel D, Lühmann D S, Sengstock K, and Weitenberg C 2016 Science 352 1091
[44]	Li T, Duca L, Reitter M, Grusdt F, Demler E, Endres M, Schleier-Smith M, Bloch I, and Schneider U 2016 Science 352 1094
[45]	Xie B Y, Wang H X, Zhang X J, Zhan P, Jiang J H, Lu M H, and Chen Y F 2021 Nat. Rev. Phys. 3 520
[46]	Po H C, Watanabe H, and Vishwanath A 2018 Phys. Rev. Lett. 121 126402
[47]	Okuma N and Sato M 2023 Annu. Rev. Condens. Matter Phys. 14 83
[48]	Herrera M A J and Bercioux D 2023 Commun. Phys. 6 42
[49]	Sun W, Yi C R, Wang B Z, Zhang W W, Sanders B C, Xu X T, Wang Z Y, Schmiedmayer J, Deng Y, Liu X J, Chen S, and Pan J W 2018 Phys. Rev. Lett. 121 250403
[50]	Wang C, Zhang P, Chen X, Yu J, and Zhai H 2017 Phys. Rev. Lett. 118 185701