Discrete Time Dirac Quantum Walk in 3+1 Dimensions

Abstract

In this paper we consider quantum walks whose evolution converges to the Dirac equation in the limit of small wave-vectors. We show exact Fast Fourier implementation of the Dirac quantum walks in one, two, and three space dimensions. The behaviour of particle states—defined as states smoothly peaked in some wave-vector eigenstate of the walk—is described by an approximated dispersive differential equation that for small wave-vectors gives the usual Dirac particle and antiparticle kinematics. The accuracy of the approximation is provided in terms of a lower bound on the fidelity between the exactly evolved state and the approximated one. The jittering of the position operator expectation value for states having both a particle and an antiparticle component is analytically derived and observed in the numerical implementations.

1. Introduction

Thinking about the discrete evolution of physical systems, the most natural example is certainly a particle moving on a lattice. A (classical) random walk is exactly the description of a particle which moves in discrete time steps and with certain probabilities from one lattice position to the neighboring lattice positions. These models have gained increasing attention, showing several applications in the fields of mathematics, physics, chemistry, computer science, natural sciences, and economics [1,2,3]. A quantum version of such a random walk—denoted quantum walk (QW)—was first introduced in [4], where the motion (right or left) of a spin-$1/2$ particle is decided by a measurement of the z-component of its spin. Subsequently, the measurement was replaced by a unitary operator on the internal space, also known as coin space [5,6], determining the evolution of the internal degree of freedom of the system. This model, known as coined quantum walk, has been proven to provide a computational speedup over classical random walks for a class of problems—such as some oracular problems, element distinctness problem, and the triangle finding problem. The Grover’s search algorithm can also be implemented as a QW [7,8,9,10,11,12,13]. The rigorous definition of QW can be found in Refs. [6,14] for the one-dimensional case, and in [5] for QW s on graphs of any dimension (see also [15] for a complete review including walks with continuous time evolution not considered in the present context).

Aside from the interest in quantum algorithms, QW s provide a fully quantum model of evolution for a system with an internal degree of freedom. As such, QW s have been considered as discrete quantum simulators for particle-physics. Interestingly, it has been proven that QW s have the capability of simulating free relativistic particle dynamics [16,17,18,19,20,21,22,23,24,25,26,27,28,29], providing—in contrast with other discretisation schemes based on finite-differences and which in general do not preserve the norm—a local unitary model underlying relativistic dynamics.

In the light of this success, in Ref. [30] the authors propose a discrete theory for quantum field dynamics based on finite dimensional quantum systems in interaction. Assuming the locality, homogeneity, and unitarity of the interaction, it follows that the systems must evolve according to a QW. Moreover, the above assumptions are very restrictive and the only QW s admissible on the cubic lattice in one, two, and three dimensions are proved to recover the usual relativistic Weyl equation in the limit of small wave-vectors. The massive case is obtained coupling two massless QW s, and in other works also the Maxwell equation [31] for Bosonic fields is proved to be compatible with an elementary QW model. Finally, the Lorentz covariance, which is broken by the discreteness of the walk, can be recovered as an approximated symmetry [32] in the relativistic limit. These results show how QW s not only provide a useful way of simulating relativistic free evolution ,but also can be considered as a fundamental approach to quantum field theory (see Refs. [33,34] for a review).

Here we consider the Dirac QW s derived in Ref. [30] and present both an analytical and a numerical study of their kinematics, recovering the characteristic traits of the usual Dirac equation.

We show “smooth-states” peaked around some wave-vector eigenstate of the QW can be considered as particle states. We present an analytical approximation of particle states evolution deriving a wave-vector-dependent differential equation for the walk evolution. Then we analyse in detail dynamical quantities such as the walk position and velocity operators and study their evolution. An intrinsic relativistic quantum processes of the Dirac field, denoted Zitterbewegung, first considered by Schrödinger [35] and corresponding to a jittering of the mean position for a relativistic particle, is recovered from the QW evolution. The theoretical existence of the quivering motion has been evidenced by numerical simulations of the Dirac equation and of quantum field theory. While Zitterbewegung oscillations cannot be directly observed by current experimental techniques for a Dirac electron since the amplitude should by very small (equal to the Compton wavelength $\hslash /mc$ with m the rest mass of the relativistic particle, namely $\approx {10}^{-12}$ m for an electron), solid state and atomic physics provide physical hardware to simulate the phenomenon [36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52], and they have recently boosted a renewed interest in the Dirac equation features.

2. Quantum Walks

A quantum walk is a local unitary evolution of a quantum system with Hilbert space $\mathcal{H}={\ell }^2\left(V\right)\otimes {\mathbb{C}}^s$, where V is a countable set and ${\mathbb{C}}^s$ is called coin space—namely, the internal degree of freedom of the walker, $s>0$ integer.

A QW on $\mathcal{H}$ is defined by assigning a mapping $\mathcal{E}:V\times V\to {\mathcal{M}}_s\left(\mathbb{C}\right)$, such that $\mathcal{E}(x,y):=U_{y,x}$, which associates to each pair of vertices a matrix, called transition matrix, acting on the coin space. Then $\psi :\mathbb{N}\to \mathcal{H}$ is a solution of the QW $(V,\mathcal{E})$ if it satisfies the following update rule for a given initial condition $\psi \left(0\right)\in \mathcal{H}$:

$$\psi (x,t+1)=\sum _{y\in V}\mathcal{E}(y,x)\psi (y,t),\forall x\in V,\forall t\in \mathbb{N},$$

(1)

where $\psi (x,t)\in {\mathbb{C}}^s$ is the QW wave-function. Since the QW evolution is unitary, the transition matrices should satisfy the following conditions for all $x,y\in V$:

$$\begin{array}{c}\hfill \sum _{z\in V}\mathcal{E}(z,x)\mathcal{E}{(z,y)}^{†}=\sum _{z\in V}\mathcal{E}{(x,z)}^{†}\mathcal{E}(y,z)={\delta }_{xy}{I}_s,\end{array}$$

(2)

where $I_s$ is the identity on the coin space ${\mathbb{C}}^s$. Such a QW carries an associated graph defined by the set of non-null transition matrices as the directed graph $\mathrm{\Gamma }=(V,E)$ with vertex set V and edge set $E:=\{(x,y)\in V\times V\mid A_{y,x}\ne 0\}$. The locality condition amounts to requiring that for every vertex $x\in V$, the cardinality of its out-neighbourhood ${\mathcal{N}}_x^{+}=\{y\in V\mid U_{y,x}\ne 0\}$, and in-neighbourhood ${\mathcal{N}}_x^{-}=\{y\in V\mid U_{x,y}\ne 0\}$, is uniformly bounded over V—namely ${\mathcal{N}}_x^{\pm }<M<+\infty$ for every $x\in V$.

In Ref. [30] it has been shown that assuming the QW homogeneous (the vertices of the graph cannot be distinguished by the walk dynamics) the graph Γ is actually a Cayley graph of a group G. Given a group G and taking $S\subseteq G$, the Cayley graph $\mathrm{\Gamma }(G,S)$ of G with connection set S is defined as the coloured directed graph $(G,S,E)$ with vertex set G, edge set $E:=\{(g,gh)\mid g\in G,h\in S\}$ and colouring given by $E\ni (g,g^{\prime })\mapsto g^{-1}{g}^{\prime }\in S$. We will assume hereafter that the connection set S is a generating set for G—which entails that the Cayley graph $\mathrm{\Gamma }(G,S)$ is unilaterally connected—and it is symmetric. Namely, $S=S^{-1}$. The walk unitary operator corresponding to the update rule of Equation (1) can be expressed in terms of the right-regular representation of G on ${\ell }^2\left(G\right)$ defined as the map $G\ni g\mapsto T_g\in \mathrm{Aut}{\ell }^2\left(G\right)$ such that $T_g|g^{\prime }\rangle =|g^{\prime }{g}^{-1}\rangle$. Assuming, by homogeneity, that we can choose the transition matrices independently of the vertex so that $U_{g,gh}\equiv U_h$ for every $g\in G$ and $h\in S$, we can write the walk operator $U\in \mathrm{Aut}\mathcal{H}$ as

$$U=\sum _{h\in S}{T}_h\otimes U_h.$$

(3)

Now the unitarity conditions on U translate into the following conditions of the transition matrices $U_h$:

$$\begin{array}{c}\hfill \sum _{h\in S}{U}_h^{†}{U}_h=\sum _{h\in S}{U}_h{U}_h^{†}=I_s,\underset{h^{-1}{h}^{\prime }=h^{\prime \prime }}{\sum _{h,h^{\prime }\in S}}{U}_h^{†}{U}_{h^{\prime }}=\underset{h^{\prime }{h}^{-1}=h^{\prime \prime }}{\sum _{h,h^{\prime }\in S}}{U}_{h^{\prime }}{U}_h^{†}=0.\end{array}$$

(4)

Fourier Representation of Abelian QW s

As pointed out by Ambainis et al. [6], there are two general ways to study the evolution of a QW. On the one hand, one can exploit the algebraic properties of the walk transition matrixes to obtain a path-sum solution, where the QW transition amplitude to a given site is expressed as a combinatorial sum over all the paths leading to that site. Regarding this approach, in Ref. [6] the authors provided a solution for the Hadamard Walk, whereas Konno derived the solution for an arbitrary coined QW [53]. Considering the application of QW s to the description of relativistic particles, also the Dirac QW in $1+1$-dimensions and the massless Dirac QW in $2+1$-dimensions have been analytically solved in position space [54,55]. On the other hand, when a QW is defined on the Cayley graph of an Abelian group, the walk dynamics can be studied in its Fourier representation, providing analytical solutions and also approximate asymptotic solutions in the long-time limit.

Let us now consider QW s defined on Cayley graphs of free Abelian groups. That is, $G\cong {\mathbb{Z}}^d$ with generating set S. Adopting the usual additive notation for the group operation on ${\mathbb{Z}}^d$, the right-regular representation of ${\mathbb{Z}}^d$ is expressed as

$$T_{\mathbf{y}}|\mathbf{x}\rangle =|\mathbf{x}-\mathbf{y}\rangle ,\forall \mathbf{x},\mathbf{y}\in {\mathbb{Z}}^d.$$

(5)

Moreover it decomposes in one-dimensional irreducible representations, as can be easily seen from the fact that the translations $T_{\mathbf{x}}$ are diagonal on the plane waves

$$|\mathbf{k}\rangle :=\frac{1}{\sqrt{\mathsf{B}}}\sum _{\mathbf{x}\in {\mathbb{Z}}^d}{e}^{-i\mathbf{k}·\mathbf{x}}|\mathbf{x}\rangle ,T_{\mathbf{x}}|\mathbf{k}\rangle =e^{-i\mathbf{k}·\mathbf{x}}|\mathbf{k}\rangle ,$$

(6)

where $\mathsf{B}$ denotes the first Brillouin zone, which depends on the specific Cayley graph $\mathrm{\Gamma }({\mathbb{Z}}^d,S)$ we are considering.

Therefore, we can write the walk operator of Equation (3) in the direct integral decomposition

$$U={\int }_{\mathsf{B}}^{\oplus }\left|\mathbf{k}\right.〉\left.〈\mathbf{k}\right|\otimes U_{\mathbf{k}}\mathbf{d}\mathbf{k},U_{\mathbf{k}}:=\sum _{\mathbf{h}\in S}{e}^{-i\mathbf{k}·\mathbf{h}}{U}_{\mathbf{h}}.$$

(7)

For each $\mathbf{k}$ we can diagonalise the matrix $U_{\mathbf{k}}$ obtaining

$$U_{\mathbf{k}}|u_r\left(\mathbf{k}\right)\rangle =e^{-i{\omega }_r\left(\mathbf{k}\right)}|u_r\left(\mathbf{k}\right)\rangle ,$$

(8)

where ${\omega }_r\left(\mathbf{k}\right)$ is the dispersion relation of the walk and $|u_r\left(\mathbf{k}\right)\rangle \in {\mathbb{C}}^s$ is the eigenvector of $U_{\mathbf{k}}$ corresponding to the eigenvalue $e^{-i{\omega }_r\left(\mathbf{k}\right)}$, with $r=1,\cdots ,s$. We notice that we have considered the representation given by the factorized orthonormal basis $|\mathbf{x}\rangle \otimes |r\rangle$ for the walk Hilbert space ${\ell }_2\left({\mathbb{Z}}^d\right)\otimes {\mathbb{C}}^s$, with ${\left\{|r\rangle \right\}}_{r=1}^s$ the canonical basis in ${\mathbb{C}}^s$.

3. The Dirac QW in One, Two, and Three Space Dimensions

In This Section, we present the Dirac QW s in one, two, and three space dimensions derived in Ref. [30]. We will see that in the limit of small wave-vectors, the Dirac walks simulate the usual Dirac equation evolution. We start from the simplest case of massless Dirac QW, also denoted Weyl QW. The massive walk will be given by coupling two Weyl QW s with the coupling parameter interpreted a posteriori as the mass of the Dirac field.

3.1. The Weyl Quantum Walk

In Ref. [30], the authors derive the unique QW s on Cayley graphs of ${\mathbb{Z}}^d$ for $d=1,2,3$ satisfying—besides locality and unitarity—the assumptions of homogeneity and discrete isotropy, and with minimal dimension s of the coin space to have non-identical evolution. As first noticed by Meyer [18], the only solution for scalar QW s on Cayley graphs of free-Abelian groups is the identical QW; in order to have non-trivial dynamics, one has to take at least $s=2$.

Let us start from Cayley graphs of ${\mathbb{Z}}^3$, the most relevant from the physical perspective. It can be proved (see [30]) that only the body-centred cubic lattice (BCC) allows one to define a QW satisfying the above assumptions. The BCC lattice is the Cayley graph $\mathrm{\Gamma }({\mathbb{Z}}^3,S_{+}\cup S_{-})$ where $S_{+}=\{{\mathbf{h}}_1,{\mathbf{h}}_2,{\mathbf{h}}_3,{\mathbf{h}}_4\}$ is the set of generators of the group and $S_{-}$ is the corresponding set of their inverses; a convenient choice for the generators is the following:

$${\mathbf{h}}_1=\left(\begin{array}{c}1\\ 1\\ 1\end{array}\right),{\mathbf{h}}_2=\left(\begin{array}{c}1\\ -1\\ -1\end{array}\right),{\mathbf{h}}_3=\left(\begin{array}{c}-1\\ 1\\ -1\end{array}\right),{\mathbf{h}}_4=\left(\begin{array}{c}-1\\ -1\\ 1\end{array}\right).$$

(9)

The first Brillouin zone $\mathsf{B}$ of the BCC lattice is defined in Cartesian coordinates as $-\pi \le k_i\pm k_j\le \pi ,i\ne j,i,j\in \{x,y,z\}$ and it is depicted in Figure 1.

The unique solutions on the BCC lattice can be summarised as:

$$\begin{array}{cccc}& U_{\mathbf{k}}=u_{\mathbf{k}}I-i\sigma ·{\tilde{\mathbf{n}}}_{\mathbf{k}},\hfill & & u_{\mathbf{k}}:=c_x{c}_y{c}_z+s_x{s}_y{s}_z,\hfill \\ & {\tilde{\mathbf{n}}}_{\mathbf{k}}:=\left(\begin{array}{c}{s}_x{c}_y{c}_z-c_x{s}_y{s}_z\\ c_x{s}_y{c}_z+s_x{c}_y{s}_z\\ c_x{c}_y{s}_z-s_x{s}_y{c}_z\end{array}\right),\hfill & & c_i:=cos{k}_i,s_i:=sin{k}_i.\hfill \end{array}$$

(10)

where σ is the vector with components given by the Pauli matrices ${\sigma }_x$, ${\sigma }_y$ and ${\sigma }_z$. The walk matrix $U_{\mathbf{k}}$ has spectrum $\{e^{-i{\omega }_{\mathbf{k}}},e^{i{\omega }_{\mathbf{k}}}\}$ with dispersion relation ${\omega }_{\mathbf{k}}=arccos{u}_{\mathbf{k}}$ and group velocity ${\mathbf{v}}_{\mathbf{k}}:={\nabla }_{\mathbf{k}}{\omega }_{\mathbf{k}}$, representing the speed of a wave-packet peaked around the central wave-vector $\mathbf{k}$.

Let us consider now $d=2$; also in this case our assumptions single out only one Cayley graph of ${\mathbb{Z}}^2$, the square lattice, involving two generators $S_{+}=\{{\mathbf{h}}_1,{\mathbf{h}}_2\}$, with ${\mathbf{h}}_1=(1,0)$ and ${\mathbf{h}}_2=(0,1)$; the first Brillouin zone $\mathsf{B}$ in this case is given by $-\pi \le k_i\le \pi ,i\in \{x,y\}$, where $k_x=k_1+k_2$ and $k_y=k_1-k_2$. The unitary matrix of the walk in Fourier representation is given by:

$$\begin{array}{cccc}\hfill & U_{\mathbf{k}}=u_{\mathbf{k}}I-i\sigma ·{\tilde{\mathbf{n}}}_{\mathbf{k}},\hfill & & u_{\mathbf{k}}:=c_x{c}_y,\hfill \\ & {\tilde{\mathbf{n}}}_{\mathbf{k}}:=\left(\begin{array}{c}{s}_x{c}_y\\ c_x{s}_y\\ s_x{s}_y\end{array}\right),\hfill & & c_i:=cos{k}_i,s_i:=sin{k}_i,\hfill \end{array}$$

(11)

with dispersion relation ${\omega }_{\mathbf{k}}=arccos{u}_{\mathbf{k}}$.

Finally, for $d=1$, the unique Cayley graph satisfying our requirements for $\mathbb{Z}$ is the lattice $\mathbb{Z}$ itself, considered as the free Abelian group on one generator $S_{+}=\left\{h\right\}$. From the unitarity conditions one gets the unique solution

$$\begin{array}{c}\hfill U_k=u_{k}I-i\sigma ·{\tilde{\mathbf{n}}}_k,u_{k}cosk,{\tilde{\mathbf{n}}}_k\left(\begin{array}{c}0\\ 0\\ sink\end{array}\right),\end{array}$$

(12)

with dispersion relation ${\omega }_k=k$.

From Equations (10) to (12), we see that the Weyl QW in dimension $d\le 3$ is of the form

$$\begin{array}{c}\hfill W_{\mathbf{k}}=u_{\mathbf{k}}I-i\sigma ·{\tilde{\mathbf{n}}}_{\mathbf{k}},\end{array}$$

(13)

for certain $u_{\mathbf{k}}$ and ${\mathbf{n}}_{\mathbf{k}}$ with dispersion relation

$$\begin{array}{c}\hfill {\omega }_{\mathbf{k}}=arccos{u}_{\mathbf{k}}.\end{array}$$

(14)

Now it is easy to show that the evolution of the walks Equations (10) to (12) obeys Weyl’s equation in the limit of small wave-vectors, and thus we call them Weyl QW s. Let us introduce the interpolating Hamiltonian $H^W\left(\mathbf{k}\right)$ defined in the wave-vector space as the matrix such that $W_{\mathbf{k}}=e^{-i{H}^W\left(\mathbf{k}\right)}$ and governing the continuous-time evolution, interpolating exactly the discrete dynamics of the walk. As one can check, the interpolating Hamiltonian is

$$\begin{array}{c}\hfill H^W\left(\mathbf{k}\right)={\omega }_{\mathbf{k}}\sigma ·{\mathbf{n}}_{\mathbf{k}},{\mathbf{n}}_{\mathbf{k}}:=\frac{1}{sin{\omega }_{\mathbf{k}}}{\tilde{\mathbf{n}}}_{\mathbf{k}},\end{array}$$

(15)

and by power expanding at the first order in $\mathbf{k}$, one has

$$\begin{array}{c}\hfill H^W\left(\mathbf{k}\right)=\sigma ·\mathbf{k}+{\mathcal{O}\left(\right|\mathbf{k}|}^2)\end{array}$$

(16)

whose first order term $\sigma ·\mathbf{k}$ coincides with the usual Weyl Hamiltonian in d dimensions, with $d=1,2,3$, once the wave-vector $\mathbf{k}$ is interpreted as the momentum.

It will be useful for the considerations of the following sections to consider the eigenvectors of the QW. Since the structure of the matrix is independent of the dimension, we will give here the general expression of the eigenvectors. Let us now rewrite the unitary matrix $W_{\mathbf{k}}$ as:

$$\begin{array}{c}\hfill W_{\mathbf{k}}=\left(\begin{array}{cc}{z}_{\mathbf{k}}& -w_{\mathbf{k}}^{*}\\ w_{\mathbf{k}}& z_{\mathbf{k}}^{*}\end{array}\right),\end{array}$$

(17)

where $z_{\mathbf{k}}$ and $w_{\mathbf{k}}$ are related to the functions in Equation (13) by the equations $Re\left(z_{\mathbf{k}}\right)=u_{\mathbf{k}}$ and ${\tilde{\mathbf{n}}}_{\mathbf{k}}=(-Im\left(w_{\mathbf{k}}\right),Re\left(w_{\mathbf{k}}\right),-Im\left(z_{\mathbf{k}}\right))$. We can then solve the eigenvalue problem

$$\begin{array}{c}\hfill W_{\mathbf{k}}|u_s^W\left(\mathbf{k}\right)\rangle =e^{-is{\omega }_{\mathbf{k}}}|u_s^W\left(\mathbf{k}\right)\rangle ,\end{array}$$

(18)

obtaining the following expression for the eigenvectors $|u_s^W\left(\mathbf{k}\right)\rangle$, with $s=\pm$,

$$\begin{array}{c}\hfill |u_s^W\left(\mathbf{k}\right)\rangle =\frac{1}{\sqrt{2}}\left(\begin{array}{c}\sqrt{1-s{v}_{\mathbf{k}}^W}\\ -s{e}^{i\phi }\sqrt{1+s{v}_{\mathbf{k}}^W}\end{array}\right),v_{\mathbf{k}}^W=\frac{Im\left(z_{\mathbf{k}}\right)}{\sqrt{1-u_{\mathbf{k}}^2}},\end{array}$$

(19)

where $\phi =Arg{w}_{\mathbf{k}}-\frac{\pi }{2}$.

3.2. The Massive Case

Now we present QW s which manifest Dirac dynamics. We consider a walk resulting from the local coupling of two Weyl QW s. One can show [30] that there is only one possible local coupling of two Weyl QW s, and that in the small wave-vector limit the resulting walks approximate the Dirac’s equation. The unique local coupling of Weyl’s QW s, modulo unitary conjugation, is of the form

$$D_{\mathbf{k}}=\left(\begin{array}{cc}n{W}_{\mathbf{k}}& im\\ im& n{W}_{\mathbf{k}}^{†}\end{array}\right),n,m\in {\mathbb{R}}^{+},n^2+m^2=1.$$

(20)

We can provide a convenient expression of the walk in terms of the gamma matrices in spinorial representation:

$$\begin{array}{c}\hfill D_{\mathbf{k}}=n{u}_{\mathbf{k}}I-in{\gamma }_0\gamma ·{\tilde{\mathbf{n}}}_{\mathbf{k}}+im{\gamma }_0,\end{array}$$

(21)

where $u_{\mathbf{k}}$ and ${\tilde{\mathbf{n}}}_{\mathbf{k}}$ are those given previously for the Weyl’s QW s. From Equation (21) we can see that the dispersion relation in this case is simply given by

$${\omega }_{\mathbf{k}}=arccos*\sqrt{1-m^2}{u}_{\mathbf{k}}.$$

(22)

In this case, the interpolating Hamiltonian $H^D\left(\mathbf{k}\right)$ has the form

$$\begin{array}{c}\hfill H^D\left(\mathbf{k}\right)=\frac{{\omega }_{\mathbf{k}}}{sin{\omega }_{\mathbf{k}}}(n{\gamma }_0\gamma ·{\tilde{\mathbf{n}}}_{\mathbf{k}}-m{\gamma }_0),\end{array}$$

(23)

and to the first order in $\mathbf{k}$ and m one obtains the usual Dirac’s Hamiltonian

$$\begin{array}{c}\hfill H^D\left(\mathbf{k}\right)={\gamma }_0\gamma ·\mathbf{k}+m{\gamma }_0+\mathcal{O}\left(m^2\right)+{\mathcal{O}\left(\right|\mathbf{k}|}^2).\end{array}$$

(24)

We notice that the Dirac QW walk not only provides the usual Dirac dispersion relation in the limit small wave-vectors, but also the correct spinorial dynamics of the Dirac equation.

It is worth noticing that in dimension $d=1$, the Dirac QW decouples into two identical $s=2$ massive QW s [22,30], written explicitly as

$$\begin{array}{c}\hfill D_k=\left(\begin{array}{cc}n{e}^{-ik}& im\\ im& n{e}^{ik}\end{array}\right),{\omega }_k=arccos*ncosk,\end{array}$$

(25)

where $n,m\in {\mathbb{R}}^{+}$, $n^2+m^2=1$. In one space dimensions, similar QW have been studied in the literature. For example in Refs. [19,28], the authors consider the relation between arbitrary coined QWs and relativistic dynamics.

For the massive QW s, the eigenvalue equation takes the form:

$$\begin{array}{c}\hfill D_{\mathbf{k}}|u_{s,p}^D\left(\mathbf{k}\right)\rangle =e^{-is{\omega }_{\mathbf{k}}}|u_{s,p}^D\left(\mathbf{k}\right)\rangle ,\end{array}$$

(26)

and the four eigenvectors $|u_{s,p}^D\left(\mathbf{k}\right)\rangle$, with $s,p=\pm$ can be written as

$$\begin{array}{c}\hfill |u_{s,p}^D\left(\mathbf{k}\right)\rangle =\frac{1}{2}\left(\begin{array}{c}\sqrt{(1-p{v}_{\mathbf{k}}^W)(1+sp{v}_{\mathbf{k}}^D)}\\ -p{e}^{i\phi }\sqrt{(1+p{v}_{\mathbf{k}}^W)(1+sp{v}_{\mathbf{k}}^D)}\\ -s\sqrt{(1-p{v}_{\mathbf{k}}^W)(1-sp{v}_{\mathbf{k}}^D)}\\ sp{e}^{i\phi }\sqrt{(1+p{v}_{\mathbf{k}}^W)(1-sp{v}_{\mathbf{k}}^D)}\end{array}\right),v_{\mathbf{k}}^D=\frac{n\sqrt{1-u_{\mathbf{k}}^2}}{\sqrt{1-n^2{u}_{\mathbf{k}}^2}},\end{array}$$

(27)

with φ, $u_{\mathbf{k}}$, and $v_{\mathbf{k}}^W$ defined for the corresponding massless QW of Equations (18) and (19).

4. Numerical Simulation of the Weyl and Dirac QW s

In order to evaluate numerically the evolution of QW s, one can adopt two different approaches. On the one hand, one can exploit the update rule in position space given by Equation (1), which is straightforward to implement numerically. This approach, however, is not very efficient if we only want to know the evolved state at some specific time t, since it would require t successive updates of the state. On the other hand, the Fourier representation of the walk allows one to directly compute the evolution at a specific time—the complexity of the computation being that of the Fourier Transform, which can be efficiently implemented via a Fast Fourier Transform (FFT) algorithm such as the Cooley–Tukey FFT algorithm [56].

Recalling the general expressions in Equations (7) and (8), the evolution of a state $|\psi \left(0\right)\rangle \in \mathcal{H}={\ell }^2\left({\mathbb{Z}}^d\right)\otimes {\mathbb{C}}^s$ is given by the subsequent application of the walk unitary $|\psi \left(t\right)\rangle =U^t|\psi \left(0\right)\rangle$. Therefore, the state at time t can be expressed in terms of its representation in Fourier space as the Fourier Transform

$$|\psi (\mathbf{x},t)\rangle =\frac{1}{{\left(2\pi \right)}^{d/2}}\sum _{r=1}^s{\int }_{\mathsf{B}}{e}^{-i\mathbf{k}·\mathbf{x}}{e}^{-i{\omega }_r\left(\mathbf{k}\right)t}{\widehat{\psi }}_r\left(\mathbf{k}\right)|u_r\left(\mathbf{k}\right)\rangle \mathbf{d}\mathbf{k},|\psi (\mathbf{x},t)\rangle \in {\mathbb{C}}^s,\mathbf{x}\in {\mathbb{Z}}^d,$$

(28)

where ${\widehat{\psi }}_r\left(\mathbf{k}\right)={\sum }_{\mathbf{x}\in G}\langle u_r\left(\mathbf{k}\right)|\psi (\mathbf{x},0)\rangle e^{i\mathbf{k}·\mathbf{x}}$ is the r-component in the eigenbasis of the walk of the discrete-time Fourier transform of $\psi (\mathbf{x},0)$. The notation of the eigenbasis refers here to that of Equation (8), where r runs over $\{1,\cdots ,s\}$ and s is the dimension of the coin. Now, the numerical data used to represent the state in Fourier space constitute a discrete sampling of it, say at frequencies $\frac{2\pi }{N_i}{k}_i$ with $k_i=-N_i/2,\cdots ,N_i/2-1$ and $N_i$ the total number of samples in dimension i. Implementing periodic boundary conditions, this amounts to taking samples in direct space over a finite region, extending the data periodically to the whole lattice.

Let us consider now the simple cubic lattice of ${\mathbb{Z}}^d$. Let us consider a restriction $f:{\mathbb{Z}}^d\to \mathbb{C}$ of $\phi \in {\ell }^2\left({\mathbb{Z}}^d\right)$ to a finite region $\mathcal{N}=\{\mathbf{m}\in {\mathbb{Z}}^d\mid 0\le m_i<N_i,i=1,\cdots ,d\}$, with $N_i\in \mathbb{N}$, such that ${f|}_{\mathcal{N}}{=\phi |}_{\mathcal{N}}$. Then, ${f|}_{\mathcal{N}}$ is periodically extended to ${\mathbb{Z}}^d$; namely $f_{\mathbf{n}+\mathbf{N}\mathbf{r}}=f_{\mathbf{n}},\forall \mathbf{n},\mathbf{r}\in {\mathbb{Z}}^d$, with periodicity matrix given by $\mathbf{N}=\mathrm{diag}(N_1,\cdots ,N_d)$. The Fourier Transform $\mathcal{F}$ of the sequence $f_{\mathbf{n}}$ coincides with the Discrete Fourier Transform (DFT) defined as

$${\widehat{f}}_{\mathbf{k}}=DFT\left(f\right)\left(\mathbf{k}\right)\frac{1}{\sqrt{N}}\sum _{\mathbf{n}\in \mathcal{N}}{f}_{\mathbf{n}}{e}^{-2\pi i{\mathbf{p}}^{⊺}{\mathbf{N}}^{-1}\mathbf{n}}{e}^{2\pi i{\mathbf{k}}^{⊺}{\mathbf{N}}^{-1}\mathbf{n}},\mathbf{k}\in \mathcal{N},$$

(29)

where $N=\mathcal{N}=\det \left(\mathbf{N}\right)$ and $p_i=*\frac{N_i}{2}$. The inversion formula is then given by:

$$f_{\mathbf{n}}={DFT}^{-1}\left(\widehat{f}\right)\left(\mathbf{n}\right)=e^{2\pi i{\mathbf{p}}^{⊺}{\mathbf{N}}^{-1}\mathbf{n}}\frac{1}{\sqrt{N}}\sum _{\mathbf{k}\in \mathcal{N}}{\widehat{f}}_{\mathbf{k}}{e}^{-2\pi i{\mathbf{k}}^{⊺}{\mathbf{N}}^{-1}\mathbf{n}}.$$

(30)

Here we have chosen the set of Fourier indices $\mathbf{k}$ so that the frequencies actually computed lie in the interval $[-\pi ,\pi ]$.

For the Dirac QW in $3+1$-dimensions, we have to consider instead the BCC lattice. One can show [57] that it is possible to reduce the DFT on the BCC lattice to two rectangular DFT s, allowing implementation of the DFT via usual rectangular FFT algorithms. We can describe the BCC lattice choosing as vertex set $G=2{\mathbb{Z}}^3\cup (2{\mathbb{Z}}^3+\mathbf{t})$, where $\mathbf{t}=(1,1,1)$. A suitable truncation of a sequence ${\phi }_{\mathbf{n}}$, $\mathbf{n}\in G$ to a function $f_{\mathbf{n}}$ defined on a finite set $B\subset G$ can be obtained choosing the fundamental region $B=2\mathcal{N}\cup (2\mathcal{N}+\mathbf{t})$ and periodically extending it to G. The original sequence $f_{\mathbf{n}}$ can be further split into two subsequences on the even and odd indices $f_{\mathbf{n}}^0=f_{2\mathbf{n}}$ and $f_{\mathbf{n}}^1=f_{2\mathbf{n}+\mathbf{t}}$, for all $\mathbf{n}\in \mathcal{N}$. As a consequence, these two sequences $f_{\mathbf{n}}^0$ and $f_{\mathbf{n}}^1$ are periodic with periodicity matrix $\mathbf{N}$: $f_{\mathbf{n}+\mathbf{N}\mathbf{r}}^j=f_{\mathbf{n}}^j$ for all $\mathbf{n}$ and $\mathbf{r}$ in ${\mathbb{Z}}^3$ and $j=0,1$. The Fourier Transform of $f_{\mathbf{n}}$ is defined as usual as

$$\begin{array}{c}\hfill {\widehat{f}}_{\mathbf{k}}=\mathcal{F}\left(f\right)\left(\mathbf{k}\right)\frac{1}{2\sqrt{B}}\sum _{\mathbf{n}\in B}{f}_{\mathbf{n}}{e}^{2\pi i{\mathbf{k}}^{⊺}{\left(2\mathbf{N}\right)}^{-1}\mathbf{n}},\forall \mathbf{k}\in \mathcal{K},\end{array}$$

(31)

where the set of Fourier indices can be chosen as $\mathcal{K}=\left\{\mathbf{k}\in {\mathbb{Z}}^3|-N_i\le k_i<N_i,i=1,2,3\right\}$. As shown in Ref. [57], one can exploit the geometry of the BCC lattice to reduce the DFT ${\widehat{f}}_{\mathbf{k}}$ with $\mathbf{k}\in \mathcal{K}$ to two functions ${\widehat{f}}_{\mathbf{k}}^0$ and ${\widehat{f}}_{\mathbf{k}}^1$ with $\mathbf{k}$ restricted now to $\mathcal{N}$. This allows for the computation of the DFT in terms of the usual rectangular DFT s:

$$\begin{array}{cc}\hfill {\widehat{f}}_{\mathbf{k}}^0& =\frac{1}{\sqrt{2}}*DFT\left(f^0\right)\left(\mathbf{k}\right)-a_{\mathbf{k}}DFT\left(f^1\right)\left(\mathbf{k}\right),\hfill \end{array}$$

(32)

$$\begin{array}{cc}\hfill {\widehat{f}}_{\mathbf{k}}^1& =\frac{1}{\sqrt{2}}*DFT\left(f^0\right)\left(\mathbf{k}\right)+a_{\mathbf{k}}DFT\left(f^1\right)\left(\mathbf{k}\right),\hfill \end{array}$$

(33)

with $\mathbf{k}\in \mathcal{N}$ and $a_{\mathbf{k}}=e^{\pi i{\mathbf{k}}^{⊺}{\mathbf{N}}^{-1}\mathbf{t}}$. Finally, from the two sequences ${\widehat{f}}_{\mathbf{k}}^0$ and ${\widehat{f}}_{\mathbf{k}}^1$, we can write the inversion formulae for $f_{\mathbf{n}}^0$ and $f_{\mathbf{n}}^1$ as:

$$\begin{array}{cc}\hfill f_{\mathbf{n}}^0& =\frac{1}{\sqrt{2}}{DFT}^{-1}({\widehat{f}}^0+{\widehat{f}}^1)\left(\mathbf{n}\right),\hfill \end{array}$$

(34)

$$\begin{array}{cc}\hfill f_{\mathbf{n}}^1& =\frac{1}{\sqrt{2}}{DFT}^{-1}{a}^{*}({\widehat{f}}^1-{\widehat{f}}^0)\left(\mathbf{n}\right).\hfill \end{array}$$

(35)

5. Kinematics of the Dirac QW

Here we study the kinematics of the Dirac QW presented in Section 3. We show that there exists a class of states whose evolution resembles the evolution of a particle with a given wave-vector. Their evolution can be described by an approximated differential equation with coefficients depending on the particle wave-vector. We observe that the positive and negative frequency eigenstates of the walk correspond to Dirac particle and antiparticle states.

Finally, we consider the position operator for the Dirac QW and find that the mean position of states having both positive and negative frequency components present the typical jittering phenomenon—denoted Zitterbewegung—of relativistic particles. The Zitterbewegung was first discovered by Schrödinger in 1930 [35], who pointed out that in the Dirac equation for free relativistic electrons, the velocity operator does not commute with the Dirac Hamiltonian. As a consequence, the evolution of the position operator shows—in addition to the classical motion proportional to the group velocity—a fast periodic oscillation with frequency $2m{c}^2$ and amplitude equal to the Compton wavelength $\hslash /mc$, with m the rest mass of the relativistic particle. This oscillating motion is due [58] to the interference of states corresponding to the positive and negative energies firstly appeared as solutions to the Dirac equation. The trembling is also shown to disappear with time [59] for a wave-packet particle state. The same phenomenology is recovered in the Dirac QW scenario that also presents solutions having positive and negative frequency eigenvalues.

5.1. Approximated Dispersive Differential Equation

We denote a quantum state of the walker a particle state if it is localized in a region of the lattice at a given instant of time and if the walk evolution preserves its localization. Accordingly, we take the following of particle state as a state that is narrow-banded in the wave-vector space.

Definition 1 (Particle-state). 

A particle-state $|\psi \rangle$ for the Dirac QW is a wave-packet smoothly peaked around some eigenvector $|u\left({\mathbf{k}}^{\prime }\right)\rangle$ of the walk. Namely, for a given ${\mathbf{k}}^{\prime }$

$$|\psi \rangle =\frac{1}{{\left(2\pi \right)}^{d/2}}{\int }_{\mathsf{B}}\mathbf{d}\mathbf{k}{g}_{{\mathbf{k}}^{\prime }}\left(\mathbf{k}\right)|\mathbf{k}\rangle |u\left(\mathbf{k}\right)\rangle ,$$

(36)

where $g_{{\mathbf{k}}^{\prime }}\in C_0^{\infty }\left[\mathsf{B}\right]$ is a smooth function satisfying the bound

$$\begin{array}{cc}& \frac{1}{{\left(2\pi \right)}^d}{\int }_{{\mathsf{B}}_{{\mathbf{k}}^{\prime }}({\sigma }_x,{\sigma }_y,{\sigma }_z)}\mathbf{d}\mathbf{k}{\left|g_{{\mathbf{k}}^{\prime }}\left(\mathbf{k}\right)\right|}^2\ge 1-\epsilon ,\epsilon >0,{\sigma }_i>0,i=x,y,z,\hfill \\ & {\mathsf{B}}_{{\mathbf{k}}^{\prime }}({\sigma }_x,{\sigma }_y,{\sigma }_z)=\{\mathbf{k}\in \mathsf{B}\mid |k_i-k_i^{\prime }|\le {\sigma }_i\}.\hfill \end{array}$$

(37)

In the next Proposition, we derive a dispersive differential equation governing the evolution for particle-states and which makes clear their particle behaviour. It will be convenient to work with the continuous time t, interpolating exactly the discrete walk evolution $U^t$. Accordingly, we consider $\mathbf{x}$, t to be real-valued continuous variable by extending the Fourier transform

$$\begin{array}{c}\hfill |\psi (\mathbf{x},t)\rangle =\frac{1}{{\left(2\pi \right)}^d}{\int }_{\mathsf{B}}\mathbf{d}\mathbf{k}{e}^{i\mathbf{k}·\mathbf{x}}{U}_{\mathbf{k}}^t|\psi (\mathbf{k},0)\rangle =\frac{1}{{\left(2\pi \right)}^d}{\int }_{\mathsf{B}}\mathbf{d}\mathbf{k}{e}^{i\mathbf{k}·\mathbf{x}}{e}^{-iH\left(\mathbf{k}\right)t}|\psi (\mathbf{k},0)\rangle ,\end{array}$$

(38)

to real $\mathbf{x}$, t. Since the walk is band-limited in momenta $\mathbf{k}\in \mathsf{B}$, then the continuous function $\psi (\mathbf{x},t)$ is completely defined by its value on the discrete points $(\mathbf{x},t)$ of the walk causal network (the sampling of a band-limited function is stated in the Nyquist-Shannon Sampling Theorem). However, all numerical results will be given only for the discrete t, namely for repeated applications of the walk unitary operator and for discrete lattice sites $\mathbf{x}$.

Proposition 1 (Dispersive differential equation). 

Consider the evolution of the Dirac QW of Section 3 on a particle-state as in Definition 1. Then for any positive integer n, the state at time t is given by

$$\begin{array}{c}\hfill |\psi (\mathbf{x},t)\rangle =e^{i({\mathbf{k}}^{\prime }·\mathbf{x}-{\omega }_{{\mathbf{k}}^{\prime }}t)}\left|\tilde{\phi }(\mathbf{x},t)\right.〉-\epsilon -\gamma {\Sigma }^{n+1}t-\mathcal{O}\left({\Sigma }^{n+3}\right)t,\end{array}$$

(39)

where $\left|\tilde{\phi }(\mathbf{x},t)\right.〉$ is solution of the following differential equation

$$\begin{array}{c}\hfill i{\partial }_t\left|\tilde{\phi }(\mathbf{x},t)\right.〉=\sum _{{\alpha }_x+{\alpha }_y+{\alpha }_z=n}\frac{{(-i)}^n{\omega }_{{\mathbf{k}}^{\prime }}^{({\alpha }_x,{\alpha }_y,{\alpha }_z)}}{{\alpha }_x!{\alpha }_y!{\alpha }_z!}\frac{{\partial }^n}{\partial x^{{\alpha }_x}\partial y^{{\alpha }_y}\partial z^{{\alpha }_z}}\left|\tilde{\phi }(\mathbf{x},t)\right.〉,\end{array}$$

(40)

$$\begin{array}{c}\hfill {\omega }_{{\mathbf{k}}^{\prime }}^{({\alpha }_x,{\alpha }_y,{\alpha }_z)}:={\left.\frac{{\partial }^n{\omega }_{\mathbf{k}}}{\partial k_x^{{\alpha }_x}\partial k_y^{{\alpha }_y}\partial k_z^{{\alpha }_z}}\right|}_{\mathbf{k}={\mathbf{k}}^{\prime }}\end{array}$$

(41)

and

$$\begin{array}{c}\hfill \gamma =(n+1)\sum _{{\alpha }_x+{\alpha }_y+{\alpha }_z=n+1}\frac{|{\omega }_{{\mathbf{k}}^{\prime }}^{({\alpha }_x,{\alpha }_y,{\alpha }_z)}|}{{\left(2\pi \right)}^d{\alpha }_x!{\alpha }_y!{\alpha }_z!}{\int }_{{\mathsf{B}}_{{\mathbf{k}}^{\prime }}({\sigma }_x,{\sigma }_y,{\sigma }_z)}\mathbf{d}\mathbf{k}{\left|g_{{\mathbf{k}}^{\prime }}\left(\mathbf{k}\right)\right|}^2,\Sigma =\max ({\sigma }_x,{\sigma }_y,{\sigma }_z).\end{array}$$

(42)

Proof. 

First we notice that at time t the particle-state in the momentum representation is simply $|\psi (\mathbf{k},t)\rangle =e^{-i{\omega }_{\mathbf{k}}t}|\psi (\mathbf{k},0)\rangle =e^{-i{\omega }_{\mathbf{k}}t}{g}_{{\mathbf{k}}^{\prime }}\left(\mathbf{k}\right)|u\left(\mathbf{k}\right)\rangle$, while in the position representation it is

$$\begin{array}{cc}\hfill & |\psi (\mathbf{x},t)\rangle :=e^{i({\mathbf{k}}^{\prime }\mathbf{x}-{\omega }_r\left({\mathbf{k}}^{\prime }\right)t)}|\phi (\mathbf{x},t)\rangle ,\hfill \end{array}$$

(43)

$$\begin{array}{cc}& |\phi (\mathbf{x},t)\rangle :=\frac{1}{{\left(2\pi \right)}^d}{\int }_{\mathsf{B}}\mathbf{d}\mathbf{k}{e}^{i(\mathbf{k}·\mathbf{x}-{\Omega }_{\mathbf{k}}t)}{g}_{{\mathbf{k}}^{\prime }}\left(\mathbf{k}\right)|u\left(\mathbf{k}\right)\rangle ,\hfill \\ \hfill & \mathbf{k}=\mathbf{k}-{\mathbf{k}}^{\prime },\hfill \\ \hfill & {\Omega }_{\mathbf{k}}={\omega }_{\mathbf{k}}-{\omega }_{{\mathbf{k}}^{\prime }}.\hfill \end{array}$$

(44)

Now we take the time derivative of $|\phi (\mathbf{x},t)\rangle$ and expand Ω vs. $\mathbf{k}$ in ${\mathbf{k}}^{\prime }$. The coefficients of the expansion can be regarded as derivatives with respect to the space coordinates and taken out of the integral (dominated derivative theorem), leading to the following dispersive differential equation:

$$\begin{array}{c}\hfill i{\partial }_t\left|\phi (\mathbf{x},t)\right.〉=\sum _{\left|\alpha \right|=1}^{\infty }\frac{{(-i)}^{\alpha }{\omega }_{{\mathbf{k}}^{\prime }}^{\left(\alpha \right)}}{\alpha !}\frac{{\partial }^{\left|\alpha \right|}}{\partial {\mathbf{x}}^{\alpha }}\left|\phi (\mathbf{x},t)\right.〉,{\omega }_{{\mathbf{k}}^{\prime }}^{\left(\alpha \right)}={\left.\frac{{\partial }^{\left|\alpha \right|}\omega \left(\mathbf{k}\right)}{\partial {\mathbf{k}}^{\alpha }}\right|}_{\mathbf{k}={\mathbf{k}}^{\prime }},\end{array}$$

(45)

where $\alpha =({\alpha }_x,{\alpha }_y,{\alpha }_z)$ is a multiindex and $\left|\alpha \right|={\alpha }_x+{\alpha }_y+{\alpha }_z$. If we truncate the above expansion at the nth order and denote by $\left|\tilde{\phi }(\mathbf{x},t)\right.〉$ the solution of the corresponding truncated differential equation, with the identification of the initial condition $\left|\tilde{\phi }(\mathbf{x},0)\right.〉=|\phi (\mathbf{x},0)\rangle$, we get the approximate state Equation (43) at time t

$$\begin{array}{c}\hfill \left|\tilde{\psi }(\mathbf{x},t)\right.〉=e^{i({\mathbf{k}}^{\prime }·\mathbf{x}-{\omega }_{{\mathbf{k}}^{\prime }}t)}\left|\tilde{\phi }(\mathbf{x},t)\right.〉.\end{array}$$

(46)

Using the definition of particle-state in Definition 1, one can compute the accuracy of the approximation Equation (46) in terms of the parameters ${\sigma }_x,{\sigma }_y,{\sigma }_z$, and ε, evaluating the overlap between the states Equations (43) and (46). That is,

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \left|\langle \tilde{\psi }\left(t\right)|\psi \left(t\right)\rangle \right|& =\left|\frac{1}{{\left(2\pi \right)}^d}{\int }_{\mathsf{B}}\mathbf{d}\mathbf{k}{e}^{-it({\sum }_{\left|\alpha \right|=n+1}\frac{{\omega }_{{\mathbf{k}}^{\prime }}^{\left(\alpha \right)}}{\alpha !}{(\mathbf{k}-{\mathbf{k}}^{\prime })}^{\alpha }+{\sum }_{\left|\alpha \right|=n+2}\mathcal{O}\left({(\mathbf{k}-{\mathbf{k}}^{\prime })}^{\alpha }\right))}{\left|g_{{\mathbf{k}}^{\prime }}\left(\mathbf{k}\right)\right|}^2\right|\hfill \\ & \ge \left|\frac{1}{{\left(2\pi \right)}^d}{\int }_{{\mathsf{B}}_{{\mathbf{k}}^{\prime }}({\sigma }_x,{\sigma }_y,{\sigma }_z)}\mathbf{d}\mathbf{k}{e}^{-it({\sum }_{\left|\alpha \right|=n+1}\frac{{\omega }_{{\mathbf{k}}^{\prime }}^{\left(\alpha \right)}}{\alpha !}{(\mathbf{k}-{\mathbf{k}}^{\prime })}^{\alpha }+{\sum }_{\left|\alpha \right|=n+2}\mathcal{O}\left({(\mathbf{k}-{\mathbf{k}}^{\prime })}^{\alpha }\right))}{\left|g_{{\mathbf{k}}^{\prime }}\left(\mathbf{k}\right)\right|}^2\right|\hfill \\ & -\left|\frac{1}{{\left(2\pi \right)}^d}{\int }_{\mathsf{B}\setminus {\mathsf{B}}_{{\mathbf{k}}^{\prime }}({\sigma }_x,{\sigma }_y,{\sigma }_z)}\mathbf{d}\mathbf{k}{e}^{-it({\sum }_{\left|\alpha \right|=n+1}\frac{{\omega }_{{\mathbf{k}}^{\prime }}^{\left(\alpha \right)}}{\alpha !}{(\mathbf{k}-{\mathbf{k}}^{\prime })}^{\alpha }+{\sum }_{\left|\alpha \right|=n+2}\mathcal{O}\left({(\mathbf{k}-{\mathbf{k}}^{\prime })}^{\alpha }\right))}{\left|g_{{\mathbf{k}}^{\prime }}\left(\mathbf{k}\right)\right|}^2\right|\hfill \\ & \ge \left|1-i(n+1)t\sum _{\left|\alpha \right|=n+1}\frac{{\omega }_{{\mathbf{k}}^{\prime }}^{\left(\alpha \right)}{\Sigma }^{n+1}}{{\left(2\pi \right)}^d\alpha !}{\int }_{{\mathsf{B}}_{{\mathbf{k}}^{\prime }}({\sigma }_x,{\sigma }_y,{\sigma }_z)}\mathbf{d}\mathbf{k}{\left|g_{{\mathbf{k}}^{\prime }}\left(\mathbf{k}\right)\right|}^2-\mathcal{O}\left({\Sigma }^{n+3}\right)t\right|-\epsilon \hfill \\ & \ge 1-\epsilon -\gamma {\Sigma }^{n+1}t-\mathcal{O}\left({\Sigma }^{n+3}\right)t\hfill \end{array}\end{array}$$

(47)

where

$$\begin{array}{c}\hfill \gamma =(n+1)\sum _{\left|\alpha \right|=n+1}\frac{|{\omega }_{{\mathbf{k}}^{\prime }}^{\left(\alpha \right)}|}{{\left(2\pi \right)}^d\alpha !}{\int }_{{\mathsf{B}}_{{\mathbf{k}}^{\prime }}({\sigma }_x,{\sigma }_y,{\sigma }_z)}\mathbf{d}\mathbf{k}{\left|g_{{\mathbf{k}}^{\prime }}\left(\mathbf{k}\right)\right|}^2,\Sigma =\max ({\sigma }_x,{\sigma }_y,{\sigma }_z).\end{array}$$

Therefore, the exact state $|\psi (\mathbf{x},t)\rangle$ at time t can be approximated by Equation (46) with the accuracy given by the overlap in Equation (47). ☐

This approximation fails to be accurate for a sufficiently large value of t. More precisely, if we require the overlap to satisfy $\left|\langle \tilde{\psi }\left(t\right)|\psi \left(t\right)\rangle \right|>1-\delta$, for some $\delta >0$, then for $t>\frac{\delta -\epsilon }{\gamma {\Sigma }^{n+1}}$ the approximated solution can deviate significantly from that of the QW. A typical application of the above proposition is the second order approximation of the state evolution. In that case, Equation (40) gives

$$\begin{array}{cc}& i{\partial }_t\left|\tilde{\phi }(\mathbf{x},t)\right.〉=\left[-i{\mathbf{v}}_{{\mathbf{k}}^{\prime }}·\nabla -\frac{1}{2}{\nabla }^T·{\mathbf{D}}_{{\mathbf{k}}^{\prime }}·\nabla \right]\left|\tilde{\phi }(\mathbf{x},t)\right.〉,\hfill \\ & {\mathbf{v}}_{\mathbf{k}}={\nabla }_{\mathbf{k}}{\omega }_{\mathbf{k}},{\mathbf{D}}_{\mathbf{k}}={\nabla }_{\mathbf{k}}{\nabla }_{\mathbf{k}}{\omega }_{\mathbf{k}},\hfill \end{array}$$

(48)

where ${\mathbf{v}}_{\mathbf{k}}$ and ${\mathbf{D}}_{\mathbf{k}}$ are respectively to drift vector and to the diffusion tensor for the particle state. Accordingly, the state will translate with group velocity given by the drift vector and its distribution in space will spread as described by the diffusion tensor.

In Figure 2, we show the numerical evolution (see Section 4) of a Gaussian particle-state in $3+1$ dimensions.

5.2. The Evolution of the QW Position Operator

Up to now we have considered only smooth-states (see Definition 1) whose walk evolution is well described by the approximate differential equation derived in Proposition 1. On the other hand, in the QW framework we are allowed to consider states very far apart from the smooth ones, and in the limit one can also consider perfectly localized states as $|\psi \rangle =|\mathbf{x}\rangle |\zeta \rangle$ with $\mathbf{x}\in {\mathbb{Z}}^d$ and $|\zeta \rangle ={\sum }_r{c}_r|r\rangle \in {\mathbb{C}}^4$, ${\sum }_r{\left|c_r\right|}^2=1$, where ${\left\{|r\rangle \right\}}_{r=1}^4$ denotes the ${\mathbb{C}}^4$ basis corresponding to the Dirac field representation in Equation (21). Since these states involve large momentum components, their evolution according to the QW dynamics will be very different from the one given by the Dirac equation. We can say that the QW determines different regimes with respect to a given reference scale at which the evolution deviates from the relativistic regime given by the Dirac equation [32,60]. However, in a QW context, the study of such states can give essential information regarding the dynamical properties of these models [5,6,14,15]. One can see in Figure 3 and Figure 4 the numerical evolution (see Section 4) of a perfectly localised state, according to the Dirac QW in $3+1$-dimensions.

The position operator X providing the representation $|\mathbf{x}\rangle$—namely, the operator such that $X|\mathbf{x}\rangle |\zeta \rangle =x|\mathbf{x}\rangle |\zeta \rangle$— is $X={\sum }_{\mathbf{x}\in {\mathbb{Z}}^d}\mathbf{x}(|\mathbf{x}\rangle \langle \mathbf{x}|\otimes I)$. Accordingly, the average position for an arbitrary one-particle state $|\psi \rangle ={\sum }_{\mathbf{x},r}{g}_r\left(\mathbf{x}\right)|\mathbf{x}\rangle |r\rangle$ is given by $\langle \psi \left|X\right|\psi \rangle$.

The definition of the mechanical momentum would need an interacting theory allowing momentum exchange between different particles. However, in Section 3, we have seen that for small $\mathbf{k}$ and m, the wave-vector $\mathbf{k}$ (namely the conjugated variable of $\mathbf{x}$ via the Fourier transform) corresponds to the Dirac particle momentum. Moreover, the momentum operator should correspond to the generator of translations over the lattice. Therefore, as conjugated momentum we take the following operator $\mathbf{P}=\frac{1}{{\left(2\pi \right)}^d}{\int }_{\mathsf{B}}\mathbf{d}\mathbf{k}\mathbf{k}(|\mathbf{k}\rangle \langle \mathbf{k}|\otimes I)$. We can now compute the commutator between $X_i$ and $P_j$, $i,j=x,y,z$. That is,

$$\begin{array}{cc}\hfill [X_i,P_j]& ={\delta }_{ij}\frac{1}{2\pi }\sum _{x_i}{\int }_{-\pi }^{\pi }\mathbf{d}{k}_j\sum _{y_i}{y}_i{k}_j|x_i\rangle \langle y_i|e^{-i{k}_j(x_i-y_i)}-{\delta }_{ij}\frac{1}{2\pi }{\int }_{-\pi }^{\pi }\mathbf{d}{k}_j\sum _{z_i,w_i}{z}_i{k}_j|w_i\rangle \langle z_i|e^{i{k}_j(z_i-w_i)}\hfill \\ & ={\delta }_{ij}\frac{1}{2\pi }\sum _{x_i,y_i}{\int }_{-\pi }^{\pi }\mathbf{d}{k}_j(x_i-y_i)k_j|x_i\rangle \langle y|e^{-i{k}_j(x_i-y_i)},\hfill \end{array}$$

(49)

where in the second equality it was possible to interchange the sum and the integral according to the Fubini Theorem. Integrating by parts we get

$$\begin{array}{c}\hfill \langle \psi |[X_i,P_j]|\psi \rangle =i\left(1-\frac{1}{2}\sum _r\left(|{\widehat{g}}_r{\left(\pi \right)|}^2+{\left|{\widehat{g}}_r(-\pi )\right|}^2\right)\right){\delta }_{ij},\end{array}$$

(50)

where $|\psi \rangle ={\sum }_{\mathbf{x},r}{g}_{\nu }\left(\mathbf{x}\right)|\mathbf{x}\rangle |r\rangle$ is a generic state and $g\left(\mathbf{k}\right)$ the discrete Fourier transform of $g\left(\mathbf{x}\right)$. We notice that Equation (50) differs from the usual canonical commutation relation by a boundary term, in agreement with the existence of perfectly localized states for the walk $|\mathbf{x}\rangle |\zeta \rangle ={\sum }_{\mathbf{y}r}{g}_r\left(\mathbf{y}\right)|\mathbf{y}\rangle |r\rangle$, $g_r\left(\mathbf{y}\right)=c_r{\delta }_{\mathbf{x}\mathbf{y}}$, for which the expectation value in Equation (50) vanishes. In the following evolution of the position expectation value, we will consider states having negligible boundary term in Equation (50).

The evolution of the position operator $\mathbf{X}\left(t\right)=U^{-t}\mathbf{X}{U}^t$ can be computed via the velocity and the acceleration operators derived by the commutator with the walk Hamiltonian

$$\begin{array}{c}\hfill \mathbf{V}\left(t\right)=i[H,\mathbf{X}\left(t\right)]={\int }_{\mathsf{B}}\mathbf{d}\mathbf{k}|\mathbf{k}\rangle \langle \mathbf{k}|\otimes \mathbf{V}\left(\mathbf{k}\right),\mathbf{A}\left(t\right)=i[H,\mathbf{V}\left(t\right)]={\int }_{\mathsf{B}}\mathbf{d}\mathbf{k}|\mathbf{k}\rangle \langle \mathbf{k}|\otimes \mathbf{A}\left(\mathbf{k}\right),\end{array}$$

(51)

From direct computation (and neglecting the boundary terms of the commutators), it follows

$$\begin{array}{cc}\hfill V_j\left(\mathbf{k}\right)& =\frac{sin{\omega }_{\mathbf{k}}-{\omega }_{\mathbf{k}}cos{\omega }_{\mathbf{k}}}{{\omega }_{\mathbf{k}}sin{\omega }_{\mathbf{k}}}H\left(\mathbf{k}\right){\left({\mathbf{v}}_{\mathbf{k}}\right)}_j+\frac{{\omega }_{\mathbf{k}}}{sin{\omega }_{\mathbf{k}}}n{\gamma }_0\gamma ·{\partial }_{k_j}{\tilde{\mathbf{n}}}_{\mathbf{k}},\hfill \end{array}$$

(52)

and

$$\begin{array}{c}\hfill A_j\left(\mathbf{k}\right)=2n\frac{{\omega }_{\mathbf{k}}^2}{{sin}^2{\omega }_{\mathbf{k}}}\left(n\sum _{\mu <\nu }{\gamma }_{\mu }{\gamma }_{\nu }{f}_{\mu \nu }^{\left(j\right)}-m\gamma ·{\partial }_{k_j}{\tilde{\mathbf{n}}}_{\mathbf{k}}\right)\end{array}$$

(53)

$$\begin{array}{c}\hfill f_{\mu \nu }^{\left(j\right)}\left({\left({\tilde{\mathbf{n}}}_{\mathbf{k}}\right)}_{\nu }{\partial }_{k_j}{\left({\tilde{\mathbf{n}}}_{\mathbf{k}}\right)}_{\mu }-{\left({\tilde{\mathbf{n}}}_{\mathbf{k}}\right)}_{\mu }{\partial }_{k_j}{\left({\tilde{\mathbf{n}}}_{\mathbf{k}}\right)}_{\nu }\right).\end{array}$$

(54)

Now we can derive the analytical expression of $\mathbf{X}\left(t\right)$ by doubly integrating the acceleration operator $\mathbf{A}\left(t\right)$, with $\mathbf{A}(\mathbf{k},t)=e^{iH\left(\mathbf{k}\right)t}\mathbf{A}\left(\mathbf{k}\right)e^{-iH\left(\mathbf{k}\right)t}$. A lengthy but simple computation (notice that $n_3{f}_{12}^{\left(j\right)}-n_2{f}_{13}^{\left(j\right)}+n_1{f}_{23}^{\left(j\right)}=0$) shows that ${[H\left(\mathbf{k}\right),\mathbf{A}\left(\mathbf{k}\right)]}_{+}=0$, which gives

$$\begin{array}{c}\hfill \mathbf{A}(\mathbf{k},t)=e^{2iH\left(\mathbf{k}\right)t}\mathbf{A}\left(\mathbf{k}\right).\end{array}$$

(55)

Therefore, integrating the first time we get

$$\begin{array}{cc}\hfill \mathbf{V}(\mathbf{k},t)& =\widehat{\mathbf{V}}\left(\mathbf{k}\right)+{\mathbf{Z}}^V(\mathbf{k},t),\hfill \\ \hfill \widehat{\mathbf{V}}\left(\mathbf{k}\right)& =\mathbf{V}\left(\mathbf{k}\right)-{\mathbf{Z}}^V(\mathbf{k},0),\hfill \\ \hfill {\mathbf{Z}}^V(\mathbf{k},t)& =\frac{1}{2i}{H}^{-1}\left(\mathbf{k}\right)\mathbf{A}(\mathbf{k},t),\hfill \end{array}$$

(56)

with $H^{-1}\left(\mathbf{k}\right)={\omega }_{\mathbf{k}}^{-2}H\left(\mathbf{k}\right)$, and integrating again one has

$$\begin{array}{cc}\hfill \mathbf{X}\left(t\right)& =\mathbf{X}\left(0\right)+\widehat{\mathbf{V}}t+{\mathbf{Z}}^X\left(t\right)-{\mathbf{Z}}^X\left(0\right),\hfill \\ \hfill {\mathbf{Z}}^X(\mathbf{k},t)& =-\frac{1}{4}{H}^{-2}\left(\mathbf{k}\right)\mathbf{A}(\mathbf{k},t),\hfill \end{array}$$

(57)

where

$$\begin{array}{cc}\hfill Z_j^X(\mathbf{k},t)& =-\frac{{\omega }_{\mathbf{k}}^3}{2sin{\omega }_{\mathbf{k}}^3}{e}^{2iH\left(\mathbf{k}\right)t}\left(n^2{\gamma }_0\gamma ·{\mathbf{w}}^{\left(j\right)}+nm\tilde{\mathbf{n}}·{\partial }_{k_j}\tilde{\mathbf{n}}+m^2{\gamma }_0\gamma ·{\partial }_{k_j}\tilde{\mathbf{n}}\right),\hfill \\ \hfill {\mathbf{w}}^{\left(j\right)}& \left(\begin{array}{c}{n}_3{f}_{13}^{\left(j\right)}+n_2{f}_{12}^{\left(j\right)}\\ -n_1{f}_{12}^{\left(j\right)}+n_3{f}_{23}^{\left(j\right)}\\ -n_1{f}_{13}^{\left(j\right)}+n_2{f}_{23}^{\left(j\right)}\end{array}\right).\hfill \end{array}$$

(58)

The operator $\widehat{\mathbf{V}}$ in Equations (56) and (57) is the classical component of the velocity operator which, in the Hamiltonian diagonal basis Equation (64), is proportional to the group velocity $\widehat{\mathbf{V}}\left(\mathbf{k}\right)\propto ({\sigma }_z\otimes I){\mathbf{v}}_{\mathbf{k}}$. In addition to the classical contribution $\widehat{\mathbf{V}}t$, we see that the position operator Equation (57) presents, as in the usual Dirac theory, a time-dependent component ${\mathbf{Z}}^X\left(t\right)$ and a constant shift term ${\mathbf{Z}}^X\left(0\right)$. Since have

$$\begin{array}{c}\hfill {}_k{\langle \mp |{\widehat{V}}_j\left(\mathbf{k}\right)|\pm \rangle }_k=0,{}_k{\langle \pm |Z_j^X\left(\mathbf{k}\right)|\pm \rangle }_k=0,\end{array}$$

(59)

the position operator $\mathbf{X}\left(t\right)$ Equation (57) mean value for the generic state $|\psi \rangle =|{\psi }_{+}\rangle +|{\psi }_{-}\rangle$ having both particle and positive and negative frequency components can always be written as

$$\begin{array}{cc}\hfill \langle \psi \left|\mathbf{X}\right(t\left)\right|\psi \rangle & ={\mathbf{x}}_{\psi }^{+}\left(t\right)+{\mathbf{x}}_{\psi }^{-}\left(t\right)+{\mathbf{x}}_{\psi }^{\mathrm{int}}\left(t\right)\hfill \end{array}$$

(60)

$$\begin{array}{cc}\hfill {\mathbf{x}}_{\psi }^{\pm }\left(t\right)& \langle {\psi }_{\pm }|\mathbf{X}\left(0\right)+\widehat{\mathbf{V}}t|{\psi }_{\pm }\rangle \hfill \end{array}$$

(61)

$$\begin{array}{cc}\hfill {\mathbf{x}}_{\psi }^{\mathrm{int}}\left(t\right)& 2Re*\langle {\psi }_{+}|\mathbf{X}\left(0\right)-{\mathbf{Z}}^X\left(0\right)+{\mathbf{Z}}^X\left(t\right)|{\psi }_{-}\rangle ,\hfill \end{array}$$

(62)

with Re denoting the real part. The first two terms ${\mathbf{x}}_{\psi }^{\pm }\left(t\right)$ simply correspond to the “classical” evolution of the particle and antiparticle components of the initial state $|\psi \rangle$, which evolve independently according to the classical component $\widehat{\mathbf{V}}$ of the velocity operator. The interference between positive and negative frequencies is responsible for the term ${\mathbf{x}}_{\psi }^{\mathrm{int}}\left(t\right)$ in Equation (62). Obviously, in case of $|\psi \rangle$ having only positive or negative component, the interference disappears. The additional term ${\mathbf{x}}_{\psi }^{\mathrm{int}}\left(t\right)$ consists of two contributions: a constant shift and a time dependent term.

Taking for example a superposition of particle and antiparticle states (see Definition 1),

$$\begin{array}{cc}& c_{+}|{\psi }_{+}\rangle +c_{-}|{\psi }_{-}\rangle ,|c_{+}{|}^2+{\left|c_{-}\right|}^2=1,\hfill \\ & |{\psi }_{\pm }\rangle =\frac{1}{{\left(2\pi \right)}^{d/2}}{\int }_{\mathsf{B}}\mathbf{d}\mathbf{k}{g}_{{\mathbf{k}}^{\prime }}(\mathbf{k},0)|\mathbf{k}\rangle |u_{\pm ,p}\left(\mathbf{k}\right)\rangle ,D_{\mathbf{k}}|u_{\pm ,p}\left(\mathbf{k}\right)\rangle =e^{\mp i{\omega }_{\mathbf{k}}}|u_{\pm ,p}\left(\mathbf{k}\right)\rangle ,\hfill \end{array}$$

(63)

where $|u_{\pm ,p}\left(\mathbf{k}\right)\rangle$ are the Dirac walk eigenvectors of Equation (27). One can show that the time dependent contribution is an oscillating term that for $t\to \infty$ goes to 0 as $1/\sqrt{t}$, and whose amplitude is bounded by $1/m$—say by the Compton wavelength $\hslash /mc$ in the usual dimensional units (see Ref [23] for the proof in one space dimension). Accordingly, ${\mathbf{x}}_{\psi }^{\mathrm{int}}\left(t\right)$ can be considered as the QW analogue of the so-called Zitterbewegung.

In Figure 5 and Figure 6, we show two numerical examples (see Section 4) of mean position evolution for the Dirac QW in one and three space dimensions, respectively. In the first case, one can also notice the time-damping of the jittering amplitudes.

Remark 1 (Newton–Wigner position operator evolution). 

As in QFT, one can define the Newton–Wigner position operator ${\mathbf{X}}_{NW}$ which does not mix states with positive and negative eigenvalues. Given the operator $W_{FW}$ providing the Foldy–Wouthuysen representation of the Dirac walk, namely the representation in which the Hamiltonian $H\left(\mathbf{k}\right)$ is diagonal

$$\begin{array}{c}\hfill W_{FW}={\int }_{\mathsf{B}}\mathbf{d}\mathbf{k}|\mathbf{k}\rangle \langle \mathbf{k}|\otimes W_{FW}\left(\mathbf{k}\right),W_{FW}\left(\mathbf{k}\right):\left\{|\nu \rangle \right\}\to \left\{{|u\rangle }_{\mathbf{k}}\right\},\end{array}$$

(64)

$$\begin{array}{c}\hfill W_{FW}^{-1}\left(\mathbf{k}\right)H\left(\mathbf{k}\right)W_{FW}\left(\mathbf{k}\right)=\mathrm{diag}({\omega }_{\mathbf{k}},{\omega }_{\mathbf{k}},-{\omega }_{\mathbf{k}},-{\omega }_{\mathbf{k}}),\end{array}$$

(65)

the Newton–Wigner rotated position operator is defined as

$$\begin{array}{c}\hfill {\mathbf{X}}_{NW}=W_{FW}^{-1}\mathbf{X}{W}_{FW}.\end{array}$$

(66)

As in the usual QFT, the Newton–Wigner position operator Equation (66) does not suffer the jittering of the mean position even for states having both a particle and an antiparticle component. Indeed, in this case, the velocity operator

$$\begin{array}{c}\hfill {\mathbf{V}}_{NW}\left(t\right)=i[H,{\mathbf{X}}_{NW}\left(t\right)],\mathbf{V}\left(k\right)=\widehat{\mathbf{V}}\left(\mathbf{k}\right),\end{array}$$

(67)

corresponds to the classical component of the velocity operator in Equation (56) and leads to a null acceleration $\mathbf{A}\left(t\right)=i[H,{\mathbf{V}}_{NW}\left(t\right)]=0$. By integrating Equation (67), we see that the time evolution of the Newton–Wigner position operator ${\mathbf{X}}_{NW}\left(t\right)$ is simply

$$\begin{array}{c}\hfill {\mathbf{X}}_{NW}\left(t\right)={\mathbf{X}}_{NW}\left(0\right)+\widehat{\mathbf{V}}t.\end{array}$$

(68)

6. Conclusions

The QW framework, say a lattice of quantum systems in local unitary interaction, appears to be very promising both from the information-theoretical perspective, in that QWs can be exploited to solve efficiently some search problems, and for the connection existing between the a discrete time quantum walk evolution and the relativistic equations of motion. In this paper, we analyse both numerically and analytically the properties of QWs on Abelian lattices up to $3+1$ dimensions. The Weyl QWs considered here are the only isotropic (all the directions on the lattice are equivalent) QWs admissible on Abelian lattices and with two-dimensional coin system. The QWs in one and two space dimensions are defined on the simple cubic lattice while the QW in $3+1$ dimensions is defined on the body-centered cubic lattice. As shown in Ref. [30], any other topology fails to accommodate a non-trivial QW (by trivial QW we mean a walk corresponding to the identical evolution or to a shift in a fixed direction). The only coupling of two Weyl QWs that preserves locality is then defined Dirac QW. Remarkably, the selected walks are compatible with a “large scale” relativistic dynamics.

The analytical results of this paper show that for particle states as defined in Definition 1, the Weyl and Dirac QW dynamics is well approximated by a dispersive differential equation whose drift and diffusion coefficients reduce to the usual Weyl and Dirac ones in the limit of small wave-vectors. The numerical results are the first simulations of QWs in $3+1$ dimensions and on the BCC lattice. The numerical results are given for the Dirac QW in $1+1$ and $3+1$ dimensions for different types of initial states. In $3+1$ dimensions we show the evolution of both particle states and perfectly localised states. In $1+1$ dimensions, the evolution of the superposition of positive and negative energy states for the Dirac QW produces (as depicted in Figure 5) the well-known Zitterbewegung effect of the relativistic electron. The appearance of this oscillating phenomenon is also shown for the Dirac QW in $3+1$ dimensions (see Figure 6).

As already mentioned, the QW framework can accommodate from a theoretical viewpoint a local discrete time unitary evolution as the microscopic description of relativistic particle dynamics. The last one is obtained as an approximation of the QW evolution for a specific class of quantum states; namely, states narrow-banded in small wave-vectors (see Definition 1). The same QW on arbitrary states (for example, localized states) shows a very different dynamical behaviour that cannot be interpreted as a particle evolution. While the approximation of Proposition 1 only works for narrow-banded states, the numerical analysis presented in the manuscript applies to arbitrary states.

Our results agree with other works (see for example Ref. [19]) in one space dimensions that studied the continuum limit of QWs, namely the lattice spacings and the time steps are sent to 0, in comparison with the Dirac or the Klein–Gordon equations. Here we do not take the same continuum limit but show that for specific input states the QW evolution recovers the relativistic one. Moreover, we do not consider only one space dimension, but also the two and three space dimensional case where the notion of spin becomes relevant.

Discrete time QWs provide a local and unitary evolution underlying the relativistic dynamics and do not start from a finite difference counterpart of the relativistic differential equations (or Hamiltonians). The main difference is in the notion of locality, since the locality of the Hamiltonian does not correspond to the locality of the unitary operator and vice versa. As a consequence the “effective” Hamiltonian corresponding to the Weyl (Dirac) QW differs from the usual Weyl (Dirac) finite difference Hamiltonian (see the $\mathrm{sinc}$ function that appears in Equations (15) and (23)). In the limit of small wave-vectors, the two Hamiltonians coincide, and both give the usual relativistic dynamics. However, for large wave-vectors they differ significantly.

The Weyl and Dirac QWs presented in this paper also provide an alternative way to discretize the usual Weyl and Dirac dynamics. The numerical results of this manuscript can be compared with other numerical approaches in the literature; see for example Ref. [61,62,63,64,65], where the authors adopt split-operator schemes to approximate the solutions of the Weyl and Dirac differential equations and recover the usual relativistic dynamics in the continuum limit.