# Dirac quantum walk on tetrahedra

Ugo Nzongani,¹,* Nathanaël Eon,² Iván Márquez-Martín,² Armando Pérez,³ Giuseppe Di Molfetta,⁴,† and Pablo Arrighi¹,‡

¹Université Paris-Saclay, INRIA, CNRS, ENS Paris-Saclay, LMF, 91190 Gif-sur-Yvette, France

²Aix-Marseille Université, Université de Toulon, CNRS, LIS, Marseille, France

³Departamento de Física Teórica and IFIC, Universidad de Valencia and CSIC, Dr. Moliner 50. 46100 Burjassot, Spain

⁴Aix-Marseille Université, Université de Toulon, CNRS, LIS, 13397 Marseille, France

Discrete-time Quantum Walks (QWs) are transportation models of single quantum particles over a lattice. Their evolution is driven through causal and local unitary operators. QWs are a powerful tool for quantum simulation of fundamental physics as some of them have a continuum limit converging to well-known physics partial differential equations, such as the Dirac or the Schrödinger equation. In this work, we show how to recover the Dirac equation in (3 + 1)-dimensions with a QW evolving in a tetrahedral space. This paves the way to simulate the Dirac equation on a curved spacetime. This also suggests an ordered scheme for propagating matter over a spin network, of interest in Loop Quantum Gravity where matter propagation has remained an open problem.

## I. INTRODUCTION

Discrete-time Quantum Walks (QWs) describe situations where a quantum particle is taking steps on a lattice conditioned on its internal state, typically a (pseudo) one half spin system. The particle dynamically explores a large Hilbert space associated with the positions of the lattice. With QWs, the transport is driven by discrete operations (shift and coin), which sets it apart from other lattice quantum simulation frameworks where transport typically rests on hopping between adjacent sites: all dynamic processes are discrete in space and time. QWs have been originally introduced by Godoy and Fujita [23, 24] as a model of coherent transport over periodic lattices. Similarly to classical random walk and Markov processes, they have applications in algorithmics, to express quantum algorithms exhibiting, in some cases, significant speedup over classical computing [1, 22, 31, 32, 37]. Most importantly for the present purpose, QWs have an immense track record for modeling and quantum simulating fundamental physics [3, 9, 13, 16, 19, 41]. It has been proved that their continuous limit, when it exists, coincides with one of the two most important equations in physics, namely the Dirac equation [10, 20, 39] and the Schrödinger equation [26]. These results have made it possible to use the QW as a framework to study the phenomenology of physical systems, such as Bloch oscillations [5, 30] or the propagation of matter in strong fields [8, 18, 34]. Moreover, they have recently proved useful as a fundamental building block for quantum simulating quantum field theories [21, 35] or studying phenomena such as weak ergodicity breaking and Hilbert space fragmentation [17, 36]. QWs have been studied in Bravais lattices in any dimensions, such as the triangular [6, 25] or cubic grid [15, 27], and for some family of regular planar graphs, such as the Apollonian networks [2]. Most of them have been proved to converge to the Dirac equation in the space-time continuous limit. For this reason,

their use is interesting in the study and simulation of Dirac materials and topological states [11, 40].

In this manuscript, we are interested in defining QWs over a tessellation of Euclidean space by means of tetrahedra. One of the main reasons why we have taken an interest in this issue is that a tetrahedron (3-simplex) is the minimal building block of 3D space, similarly to the role played by a triangular tessellation in a 2D space. The tessellation that we consider is also the simplest way to fill Euclidean space with tetrahedra [38]. As it was shown in [6, 25], the Dirac equation can be simulated as the continuum limit of a QW defined on a triangular or honeycomb lattice. This result was further generalized to a curved spacetime [7], where the continuum limit also describes the corresponding Dirac equation, by exploiting the duality between changes in geometry, and changes in local unitaries that incorporate the deformation of the metric due to the lack of flatness in spacetime. A crucial result from this work was that the square lattice fails in reproducing this duality for arbitrary deformations, in contrast to the triangular and honeycomb tessellations. From a more practical perspective, the latter implementations of a QW on a curved geometry need fewer steps than QWs defined on a rectangular grid [4]. All these qualities make the QW based on tetrahedra a very promising candidate to investigate the above results for the 3D space.

Finally, another reason for the tetrahedra covering of space is that it is dual to the 4-valent node, the elementary block of a spin network [12]. In Loop Quantum Gravity, a spin network is a basis state of the Hilbert space of the quantized gravitational field on a 3-dimensional hypersurface [33]. Modelling matter propagation over a spin network is still an open problem [14]. In this paper we address this problem by introducing a toy model where a quantum walker propagates on a regular spin network. We prove that we recover, in the continuous space-time limit, the Dirac equation in (3+1)-dimensions.

The organization of the manuscript is the following: in Sec. II we introduce the formal mathematical model of the QW over tetrahedra. In Sec. III we briefly recall the basics of the Dirac equation and its relation to the Weyl equation. In Sec. IV we show how to simulate the Dirac equation in 3 + 1 dimensions on a tetrahedral space. We also discuss

<footnote>* ugo.nzongani@universite-paris-saclay.fr</footnote>

<footnote>† giuseppe.dimolfetta@lis-lab.fr</footnote>

<footnote>‡ pablo.arrighi@universite-paris-saclay.fr</footnote>
how this provides a scheme whereby matter propagates over a spin network. We summarize our main results in Sec. V. Appendix B provides a robust generalisation to surfaces with boundaries, but only at the cost of doubling the number of degrees of freedom. We work in units for which $\hbar = c = 1$.

## II. TETRAHEDRAL QUANTUM WALK

Among all known ways of filling Euclidean space with tetrahedra [38], the simplest (albeit non-trivial) way is to dissect a cube into six 3-orthoschemes, i.e. tetrahedra where all four faces are right triangles. All of them have identical size. However, three of them are left-handed (LH) and the three others are right-handed (RH) (i.e. mirror to each other, since we need one LH and one RH to form a cube face) as shown in Fig. 1. At time step $t$, each tetrahedron $k$ hosts a $\mathbb{C}^4$ vector:

$$ \phi(t, k) = (\phi(t, k, i))_{i=0 \dots 3}^T, \quad (1) $$

with each complex amplitude $(\phi(t, k, i))_{i=0 \dots 3}^T$ located on a different facet of the tetrahedron, as represented in Fig. 1 by coloured dots: blue, red, cyan and magenta corresponding to facets $i = 0, 1, 2, 3$ respectively. Assuming there are no boundaries, each facet is shared by two tetrahedra, a $\mathbb{C}^2$ vector on each shared facet can be defined. In the following we focus on the regular lattice filling up 3D space with such cubes, unrotated<fnref n="1" />, as shown in Fig. 4a. This allows us to adopt the convenient convention that a LH (resp. RH) tetrahedron be identified in space by means of the coordinates $k = (k_x, k_y, k_z)^T$ of its blue $i = 0$ (resp. red $i = 1$) component with $k \in \varepsilon\mathbb{Z}^3$, with $\varepsilon$ the characteristic lattice spacing which will be defined below. Notice that, as shown in Fig. 2, a cube composed of 6 tetrahedra can be further dissected into 8 smaller and identical cubes. The complex amplitudes are then all located in the middle of the small cubes and at the center of their faces. Let $\varepsilon$ be the size of such small cubes. Thus, in a given tetrahedron the distance between the blue and red components is $\varepsilon$, and the one separating a blue (resp. red) from a cyan (resp. magenta) is $\varepsilon/2$.

Dirac QWs over the cubic lattice are of course well-known, and our strategy here will be to fall back on our feet and recover such a known scheme, in spite of the fact that the basic constituents of the lattice are the tetrahedra that make up the cubes and not the cubes themselves. This further subdivision does introduce a number of difficulties in terms of laying out the spinors and orchestrating their propagation, but more fundamentally also in terms of locality: each local unitary gate that makes up the QW needs be local to an individual tetrahedron, or to the joint facet of two glued tetrahedra. The Tetrahedral QW proposed aims to achieve a massless Dirac QW first. Later in the paper we add the mass term to obtain

<figure><img src="image_1.png" /><figcaption>FIG. 1: (Top) Left-handed and right-handed tetrahedra. (Bottom) Dissection of a cube into six tetrahedra. Coloured dots blue, red, cyan and magenta correspond to the amplitudes lying on the facets $i = 0, 1, 2, 3$ respectively.</figcaption></figure>

a Dirac QW, again in a way that respects the refined locality constraints of the tetrahedral tesselation.

Across the paper, we say that a quantum operation is causal when it is allowed to push information from one tetrahedra to its neighbour, and so on, forming a chain. We say that it is localized if it decomposes into local operations. We distinguish, however, three types of local operations. Those acting solely on a tetrahedron will be referred to as strictly local. Those acting solely on a facet will be referred to as facet-local. Those allowed to access both the internal amplitudes of a tetrahedron and some amplitudes across its facets will be referred to as weakly local.

### A. Shift operator

The shift operator $S$ decomposes in two steps defined by $S_B$ and $S_G$. $S_B$ is causal: for the LH (resp. RH) (i) the blue and cyan (resp. red and magenta) components go into the red and magenta (resp. blue and cyan) facets. At the same time,

<footnote n="1">¹ Thus, cube-cube glueings are RH-RH and LH-LH glueings. Our attempts to systematically enforce RH-LH alternation do not allow us to implement an $x \to y \to z$ sequence of displacements.</footnote>
<figure><img src="image_2.png" /><figcaption>FIG. 2: Dissection of a cube with 8 smaller cubes.</figcaption></figure>

the red (resp. blue) moves to the cyan (resp. magenta) that lies across the magenta (resp. cyan) facet. Synchronously, the magenta (resp. cyan) moves to the blue (resp. red) of its neighboring tetrahedron. $S_G$ is localized, i.e. it is implemented by strictly local operations on each tetrahedron: (ii) the blue and cyan (resp. red and magenta) components are swapped. The first and second steps correspond, respectively, to the black and grey arrows shown in Fig. 3. Let $n(k, i)$ be the neighbor of tetrahedron $k$ on facet $i$. The shifting dynamics therefore reads:

$$ \phi(t + \varepsilon, k) = (S_G S_B \phi)(t, k), \quad (2) $$

with

$$ (S_B \phi)(t, k) = \begin{cases} \begin{pmatrix} \phi(t, n(k, 2), 3) \\ \phi(t, k, 0) \\ \phi(t, n(k, 2), 1) \\ \phi(t, k, 2) \\ \phi(t, k, 1) \\ \phi(t, n(k, 3), 2) \\ \phi(t, k, 3) \\ \phi(t, n(k, 3), 0) \end{pmatrix} & \text{if } k \text{ is LH} \\ \begin{pmatrix} \phi(t, n(k, 2), 3) \\ \phi(t, k, 0) \\ \phi(t, n(k, 2), 1) \\ \phi(t, k, 2) \\ \phi(t, k, 1) \\ \phi(t, n(k, 3), 2) \\ \phi(t, k, 3) \\ \phi(t, n(k, 3), 0) \end{pmatrix} & \text{otherwise} \end{cases}, \quad (3) $$

$$ (S_G \phi)(t, k) = \begin{cases} \begin{pmatrix} \phi(t, k, 2) \\ \phi(t, k, 1) \\ \phi(t, k, 0) \\ \phi(t, k, 3) \\ \phi(t, k, 0) \\ \phi(t, k, 3) \\ \phi(t, k, 2) \\ \phi(t, k, 1) \end{pmatrix} & \text{if } k \text{ is LH} \\ \begin{pmatrix} \phi(t, n(k, 2), 1) \\ \phi(t, k, 0) \\ \phi(t, n(k, 2), 3) \\ \phi(t, k, 2) \\ \phi(t, k, 1) \\ \phi(t, n(k, 3), 0) \\ \phi(t, k, 3) \\ \phi(t, n(k, 3), 2) \end{pmatrix} & \text{otherwise} \end{cases}. \quad (4) $$

As we show on Fig. 4b, the global effect of this dynamics corresponds to the action of a single operator $S$:

$$ \phi(t + \varepsilon, k) = (S\phi)(t, k), \quad (4) $$

<figure><img src="image_3.png" /><figcaption>FIG. 3: Shifting dynamics of the tetrahedral QW working in two steps. First the black arrows act causally on the tetrahedra, with the gradient effect illustrating the need to chain them up. Second the grey arrows (thin ones) that act in a strictly localized manner. Notice that only the coin step shuffles the amplitudes between LH and RH tetrahedra.</figcaption></figure>

with

$$ (S\phi)(t, k) = \begin{cases} \begin{pmatrix} \phi(t, n(k, 2), 1) \\ \phi(t, k, 0) \\ \phi(t, n(k, 2), 3) \\ \phi(t, k, 2) \\ \phi(t, k, 1) \\ \phi(t, n(k, 3), 0) \\ \phi(t, k, 3) \\ \phi(t, n(k, 3), 2) \end{pmatrix} & \text{if } k \text{ is LH} \\ \begin{pmatrix} \phi(t, n(k, 2), 1) \\ \phi(t, k, 0) \\ \phi(t, n(k, 2), 3) \\ \phi(t, k, 2) \\ \phi(t, k, 1) \\ \phi(t, n(k, 3), 0) \\ \phi(t, k, 3) \\ \phi(t, n(k, 3), 2) \end{pmatrix} & \text{otherwise} \end{cases}. \quad (5) $$

Notice that the blue and red components are moving along the cubic grid of Fig. 2 avoiding the forbidden positions. By inspecting at Fig. 4b, we see that the cyan and magenta components are following them closely behind, as if the cyan (resp. magenta) component were linked to the blue (resp. red).
<figure><img src="image_4.png" /><figcaption>(a) The upper left (resp. right) figure show the global dynamics of the first (resp. second) walker. On the lower figures, we show the same dynamics but replacing the tetrahedra by small cubes. The origin coordinate k corresponds to the black dot located at the bottom right of the small cubes.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIG. 4: Global dynamics of the tetrahedral QW and trajectories of the first walker.</figcaption></figure>

These two walkers are independent of each other. They cross each other but do not interact. In other words, as proved in Appendix A, we have here two independent walkers evolving according to the same dynamics, the first one on the blue and red components, and the second one on the cyan and magenta components.

Let us focus on the first walker. The two-dimensional wave function living on the shared facet of neighboring tetrahedra are defined as follows:

$$ \psi(t, k) = \begin{pmatrix} \psi^{\uparrow}(t, k) \\ \psi^{\downarrow}(t, k) \end{pmatrix}, \qquad (6) $$

where the up (resp. down) component lies in the RH (resp.

LH) tetrahedra. Let us focus on what happens at coordinate k for instance:

$$ \psi(t + \varepsilon, k) = \begin{pmatrix} \psi^{\uparrow}(t, k + \varepsilon u_z) \\ \psi^{\downarrow}(t, k - \varepsilon u_z) \end{pmatrix}, \qquad (7) $$

where $u_z$ is the Bloch vector of $\sigma_z$ i.e. the unit vector along the z-axis. It is clear from Eq. (7) that the up (resp. down) component is going down (resp. up) in the 3-coordinate which for Fig. 4a is vertical. Along one axis, say z, we see that the above finite difference equations converge in the continuous limit $\varepsilon \to 0$ to the transport equation:

$$ (\partial_t \psi)(t, k) = (\sigma_z \partial_z \psi)(t, k). \qquad (8) $$
Furthermore, as the axis shown in the center of Fig. 4a are the one we use and by looking at the direction of the black arrows in Fig. 4b, we see that only black arrows that goes up (resp. down) belongs to an LH (resp. RH) tetrahedron. Thus, the upper (resp. lower) component of $\psi(t, k)$ is located in the RH (resp. LH) tetrahedra.

The direction in which the components moves depends on the coordinates, (again by looking at Fig. 4a):

$$u(k) = \begin{cases} u_x & \text{if } k = \varepsilon(\text{even}, \text{odd}, \text{odd})^T \text{ or } \varepsilon(\text{odd}, \text{even}, \text{even})^T \\ u_y & \text{if } k = \varepsilon(\text{even}, \text{odd}, \text{even})^T \text{ or } \varepsilon(\text{odd}, \text{even}, \text{odd})^T \\ u_z & \text{if } k = \varepsilon(\text{even}, \text{even}, \text{odd})^T \text{ or } \varepsilon(\text{odd}, \text{odd}, \text{even})^T \\ \emptyset & \text{if } k = \varepsilon(\text{even}, \text{even}, \text{even})^T \text{ or } \varepsilon(\text{odd}, \text{odd}, \text{odd})^T \end{cases} \quad (9)$$

where $u_\mu$ is the Bloch vector of $\sigma_\mu$ e.g. $u_x = (1, 0, 0)^T$. Inspecting Fig. 4b we see that the z-axis movers become x-axis movers. In general, $\mu$-axis movers become $(\mu + 1)$-axis movers, as made formal by $u_{\mu+1} = Pu_\mu$ with:

$$P = \begin{pmatrix} 0 & 0 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \end{pmatrix}. \quad (10)$$

The shift evolution of the walker then reads:

$$\begin{aligned} \psi(t + \varepsilon, k) &= \hat{S} \begin{pmatrix} \psi^{\uparrow}(t, k) \\ \psi^{\downarrow}(t, k) \end{pmatrix} \\ &= \begin{pmatrix} \psi^{\uparrow}(t, k + \varepsilon u(k)) \\ \psi^{\downarrow}(t, k - \varepsilon u(k)) \end{pmatrix} \quad (11) \\ &= \begin{pmatrix} \tau_{u(k), \varepsilon} \psi^{\uparrow}(t, k) \\ \tau_{u(k), -\varepsilon} \psi^{\downarrow}(t, k) \end{pmatrix}, \end{aligned}$$

where $\hat{S}$ denotes the effect of $S$ but on the two-dimensional $\psi$'s, and $\tau_{\mu, \varepsilon}$ is the translation by $\varepsilon$ along the $\mu$-axis. This looks like it does implement the correct usual partial shift $T_{u, \varepsilon}$ in the $u(k)$ direction of a cubic QW. Let us look at this in more detail. On a cubic lattice, a partial shift in some direction $u$ is written as

$$T_{u, \varepsilon} \begin{pmatrix} \psi^{\uparrow}(t, k) \\ \psi^{\downarrow}(t, k) \end{pmatrix} = \begin{pmatrix} \psi^{\uparrow}(t, k + \varepsilon u) \\ \psi^{\downarrow}(t, k - \varepsilon u) \end{pmatrix} = \begin{pmatrix} \tau_{u, \varepsilon} \psi^{\uparrow}(t, k) \\ \tau_{u, -\varepsilon} \psi^{\downarrow}(t, k) \end{pmatrix}. \quad (12)$$

### B. Coin operator

Since the walker changes direction along the dynamics, and since our aim is to mimic a transport equation such as $\partial_t \psi = \sigma_x \partial_x \psi + \sigma_y \partial_y \psi + \sigma_z \partial_z \psi$, we have to introduce basis changes to map $\sigma_{\mu-1}$ into $\sigma_\mu$. We notice that $R_{\sigma_z}(\theta) R_{\sigma_x}(\theta)$ with $\theta = \pi/2$ maps the Bloch vector of $\sigma_{\mu-1}$ into that of $\sigma_\mu$ with $\mu = 1, 2, 3$, i.e:

$$\sigma_{\mu} = R_{\sigma_z}(\theta) R_{\sigma_x}(\theta) \sigma_{\mu-1} R_{\sigma_x}^{\dagger}(\theta) R_{\sigma_z}^{\dagger}(\theta) \quad (13)$$

since $(R_{\sigma_z}(\theta) R_{\sigma_x}(\theta))^3 = -I$ and we want it to be the identity, we define $C = e^{i\frac{\pi}{3}} R_{\sigma_z}(\theta) R_{\sigma_x}(\theta)$. This $C$ operation is applied on each facet at once, it is therefore facet-localized.

<figure><img src="image_6.png" /><figcaption>FIG. 5: Independent sub-walks of a 2-dimensional walk. Forget about grey arrows for now. As the amplitudes that were at blue positions follow the black arrows, they now find themselves at red positions, without ever interacting with those which were red initially, who have now left to blue positions.</figcaption></figure>

### C. Complete dynamics

It is important to observe that, within the first quantum walker whose propagation is described by Fig. 4b, there exist in fact three sub-walkers that are “independent”, in the sense that they avoid each other and never interact. This is easier to see in 2 spatial dimensions, see Fig. 5, where the equivalent of $\hat{S}$ features two independent sub-walkers. If we remove one of the two sub-walkers, we see that the one remaining follows operator $T_x$, then $T_y$, etc. This is why $C\hat{S}$ is a suitable implementation of $CT_yCT_x$.

## III. THE DIRAC EQUATION

Let us briefly recall some properties of the Dirac equation [29], in order to better understand the connection with the QW defined above. The Dirac equation describes the motion of a relativistic spin 1/2 particle with mass $m$. It can be written as

$$(i\gamma^\mu \partial_\mu - m)\psi = 0, \quad (14)$$

where $\gamma^\mu$ ($\mu = 0, 1, 2, 3$) are the gamma matrices. This equation can be decomposed into a pair of coupled Weyl-like (i.e., with and additional mass term) equations [28] acting on the first two and last two indices of the original four-component spinor, i.e. by writing

$$\psi = \begin{pmatrix} \psi_R \\ \psi_L \end{pmatrix}, \quad (15)$$

where $\psi_L$ and $\psi_R$ are each two-component Weyl spinors. This can be achieved by using the so-called chiral representation of the gamma matrices [29]. Following this procedure, Eq. (14) becomes equivalent to a pair of coupled Weyl equations

$$i\sigma^\mu \partial_\mu \psi_R = m\psi_L \quad (16)$$
$$i\bar{\sigma}^{\mu}\partial_{\mu}\psi_{L} = m\psi_{R}. \qquad (17)$$

In the above equations, $\sigma^\mu$ has components $(I_2, \sigma^i)$ and $\bar{\sigma}^\mu$ has components $(I_2, -\sigma^i)$, with $I_2$ the two-dimensional identity matrix. In the massless $m=0$ case both equations become two true decoupled Weyl equations. The coupling between both equations (16) and (17) is associated to the action of the $\gamma^0$ matrix, which “flips” both Weyl spinors.

In the following sections we will make use of the above ideas to construct a QW that describes the massive Dirac equation. First we construct a QW for the massless Dirac equation. Then, the necessary flipping operation will be introduced to include the mass term.

# IV. DIRAC QUANTUM WALK

## A. Massless Dirac QW

After multiplication by $\gamma^0$, the massless Dirac equation takes the form

$$i\partial_0 \Psi = \hat{D}\Psi \quad \text{with} \\ \hat{D} = i \sum_{j=1}^{3} \alpha^j \partial_j, \qquad (18)$$

where $\alpha^j = \gamma^0 \gamma^j$. To obtain the massless Dirac equation, we therefore need a four-dimensional spinor and a new coin acting on these four amplitudes. Let $\Psi(t, k) = (\Psi_i(t, k))_{i=0,...,3}$ be such a spinor. Its components are the coloured dots located on the small cubes of Fig. 4b. Therefore, such a spinor contains components of two tetrahedra, LH and RH. Note that this spinor contains either blue and cyan, or red and magenta components, depending on the parity of the position. For those containing blue and cyan (resp. red and magenta) components; the first component ($i=0$) is the cyan (resp. magenta) dot of the RH tetrahedron, the second ($i=1$) is the blue (resp. red) of the RH, the third ($i=2$) is the blue (resp. red) of the LH, and the last one ($i=3$) is the cyan (resp. magenta) of the LH tetrahedron. It turns out that the chiral representation keeps the strict locality, thus:

$$\alpha^j = \sigma_3 \otimes \sigma_j. \qquad (19)$$

The coin operator $\hat{C}$ must verify $\alpha^\mu = \hat{C}\alpha^{\mu-1}\hat{C}^\dagger$. That way, $\hat{C} = \sigma_0 \otimes C$ is the coin operator. We keep identical the first (black) and second (grey) steps of the shift operation, only

reexpressing their effects in terms of these spinors:

$$(T_{B,u(k),\epsilon}\Psi)(t,k) = \begin{cases} \begin{pmatrix} \Psi_0(t, k + \epsilon u(k)) \\ \Psi_1(t, k + \epsilon u(k)) \\ \Psi_3(t, k - \epsilon u(k)) \\ \Psi_2(t, k - \epsilon u(k)) \end{pmatrix} & \text{if } k \text{ is LH} \\ \begin{pmatrix} \Psi_1(t, k + \epsilon u(k)) \\ \Psi_0(t, k + \epsilon u(k)) \\ \Psi_2(t, k - \epsilon u(k)) \\ \Psi_3(t, k - \epsilon u(k)) \end{pmatrix} & \text{otherwise} \end{cases}, \\ \\ (T_{G,u(k),\epsilon}\Psi)(t,k) = \begin{cases} \begin{pmatrix} \Psi_0(t, k) \\ \Psi_1(t, k) \\ \Psi_3(t, k) \\ \Psi_2(t, k) \end{pmatrix} & \text{if } k \text{ is LH} \\ \begin{pmatrix} \Psi_1(t, k) \\ \Psi_0(t, k) \\ \Psi_2(t, k) \\ \Psi_3(t, k) \end{pmatrix} & \text{otherwise} \end{cases}. \qquad (20)$$

Hence,

$$(T_{G,b(k),\epsilon} T_{B,u(k),\epsilon} \Psi)(t,k) = \begin{pmatrix} \Psi_0(t, k + \epsilon u(k)) \\ \Psi_1(t, k + \epsilon u(k)) \\ \Psi_2(t, k - \epsilon u(k)) \\ \Psi_3(t, k - \epsilon u(k)) \end{pmatrix}, \qquad (21)$$

performs a partial shift in the $u(k)$ direction. For the $z$ axis, this approximates as:

$$(\mathbb{I} + \epsilon(\sigma_3 \otimes \sigma_3)\partial_3) \Psi = T_{G,3,\epsilon} T_{B,3,\epsilon} \Psi + O(\epsilon^2). \qquad (22)$$

In the same way, along the $y$-axis,

$$(\mathbb{I} + \epsilon(\sigma_3 \otimes \sigma_2)\partial_2) \Psi = \hat{C}^\dagger T_{G,2,\epsilon} T_{B,2,\epsilon} \hat{C} \Psi + O(\epsilon^2). \qquad (23)$$

Lastly, for the $x$-axis we obtain:

$$(\mathbb{I} + \epsilon(\sigma_3 \otimes \sigma_1)\partial_1) \Psi = \hat{C} T_{G,1,\epsilon} T_{B,1,\epsilon} \hat{C}^\dagger \Psi + O(\epsilon^2). \qquad (24)$$

The time evolution of the massless Dirac QW thus reads:

$$\Psi(t + \epsilon, k) = \hat{W}_{\epsilon} \Psi(t, k), \qquad (25)$$

with:

$$\hat{W}_{\epsilon} = \hat{C}^{\dagger} T_{G,2,\epsilon} T_{B,2,\epsilon} \hat{C}^{\dagger} T_{G,1,\epsilon} T_{B,1,\epsilon} \hat{C}^{\dagger} T_{G,3,\epsilon} T_{B,3,\epsilon}. \qquad (26)$$
<figure><img src="image_7.png" /><figcaption>FIG. 6: Displacement graph of the tetrahedral Dirac QW.</figcaption></figure>

Using Eq. (B15) and Eq. (26) we get:

$$ \Psi(t + \varepsilon, k) = \prod_{j=1}^{3} \hat{C}^{\dagger}(\mathbb{I} + \varepsilon(\sigma_3 \otimes \sigma_3)\partial_j)\Psi(t, k) + O(\varepsilon^2) $$

$$ \Psi(t + \varepsilon, k) = \prod_{j=1}^{3} (\mathbb{I} + \varepsilon(\sigma_3 \otimes \sigma_j)\partial_j)\Psi(t, k) + O(\varepsilon^2) $$

$$ (\mathbb{I} + \varepsilon\partial_0)\Psi(t, k) = \left(\mathbb{I} + \sum_{j=1}^{3} \varepsilon(\sigma_3 \otimes \sigma_j)\partial_j\right) \Psi(t, k) + O(\varepsilon^2) $$

$$ \partial_0 \Psi(t, k) = \sum_{j=1}^{3} (\sigma_3 \otimes \sigma_j) \partial_j \Psi(t, k) + O(\varepsilon^2). \qquad (27) $$

Therefore, in the continuum limit we approximate the massless Dirac equation in (3+1)-dimensions.

## B. Dirac with mass

The QW we introduced only portrays the movement of massless particles. To prevent them from travelling at the speed of light we introduce an additional new coin acting as the mass term. The massive Dirac equation now takes the form:

$$ i\partial_0\Psi = \mathcal{D}\Psi \quad \text{with} \quad \mathcal{D} = m\gamma^0 + i\sum_{j=1}^{3} \alpha^j \partial_j, \qquad (28) $$

where $m$ is the mass of the particle. The $\alpha^j$ matrices are that of Eq. (19) and $\gamma^0$ is defined as:

$$ \gamma^0 = \sigma_1 \otimes \sigma_0. \qquad (29) $$

Let $M = e^{-i\varepsilon m \gamma^0}$ be the mass coin. According to our choice of representation in Eq. (19) and Eq. (29), the coin $\hat{C}$ makes the first and second (resp. third and fourth) components of $\Psi$ interact with each others, meaning that it is a strictly localized operator. On the other hand, the mass coin $M$ makes the first and third (resp. second and fourth) components interact with

each others, which means that it is weakly localized i.e. it can be decomposed into weakly local operators, each of which is allowed to access an amplitude from across one of the facets of the tetrahedron.

Moreover, we notice that:

$$ \begin{aligned} \hat{C}M\hat{C}^{\dagger} &= M \\ M\hat{C}M^{\dagger} &= \hat{C}, \end{aligned} \qquad (30) $$

making it is possible to successively apply the coins without disturbing the walk. The time evolution of the massive Dirac tetrahedral QW is defined as:

$$ \Psi(t + \varepsilon, k) = (W_{\varepsilon}\Psi)(t, k), \qquad (31) $$

where

$$ \begin{aligned} W_{\varepsilon} &= \hat{C}^{\dagger} \mathcal{T}_{G,2,\varepsilon} M \mathcal{T}_{B,2,\varepsilon} \hat{C}^{\dagger} \mathcal{T}_{G,1,\varepsilon} M \mathcal{T}_{B,1,\varepsilon} \hat{C}^{\dagger} \mathcal{T}_{G,3,\varepsilon} M \mathcal{T}_{B,3,\varepsilon}, \\ &= \hat{C}_{G,2,\varepsilon} \mathcal{T}_{B,2,\varepsilon} \hat{C}_{G,1,\varepsilon} \mathcal{T}_{B,1,\varepsilon} \hat{C}_{G,3,\varepsilon} \mathcal{T}_{B,3,\varepsilon} \end{aligned} \qquad (32) $$

with $\hat{C}_{G,\mu,\varepsilon} = \hat{C}^{\dagger} \mathcal{T}_{G,\mu,\varepsilon} M$.

In the end, the tetrahedral Dirac QW is just the repetition of two steps, one being causal and the other one localized. More precisely, according to Eq. (32) it consists in: (i) a causal shift $\mathcal{T}_{B,\mu,\varepsilon}$ and (ii) applying the mass coin $M$, performing a shift $\mathcal{T}_{W,\mu,\varepsilon}$ and applying a basis change coin $\hat{C}^{\dagger}$, all three localized.

In the continuum limit we do obtain the massive Dirac equation since:

$$ \begin{aligned} \Psi(t + \varepsilon, k) &= (\mathbb{I} - i\varepsilon D)\Psi(t, k) + O(\varepsilon^2) \\ &= \prod_{j=1}^{3} \hat{C}^{\dagger}(\mathbb{I} + \varepsilon(\sigma_3 \otimes \sigma_3)\partial_j)M\Psi(t, k) + O(\varepsilon^2) \\ &= \prod_{j=1}^{3} (\mathbb{I} + \varepsilon(\sigma_3 \otimes \sigma_j)\partial_j)M\Psi(t, k) + O(\varepsilon^2). \end{aligned} \qquad (33) $$

Moreover, as $M = \mathbb{I} - i\varepsilon m \gamma^0 + O(\varepsilon^2)$, we have:

$$ \begin{aligned} (\mathbb{I} + \varepsilon\partial_0)\Psi(t, k) &= \left(\mathbb{I} - 3im\varepsilon\gamma^0 + \sum_{j=1}^{3} \varepsilon(\sigma_3 \otimes \sigma_j)\partial_j\right) \Psi(t, k) \\ &\quad + O(\varepsilon^2) \\ \partial_0 \Psi(t, k) &= \left(-3im\gamma^0 + \sum_{j=1}^{3} (\sigma_3 \otimes \sigma_j)\partial_j\right) \Psi(t, k) \\ &\quad + O(\varepsilon^2). \end{aligned} \qquad (34) $$

Let us notice that it would be possible to modify the QW operator defined in Eq. (32) so that the mass term appears only once, by redefining

$$ W_{\varepsilon} = M \hat{C}^{\dagger} \mathcal{T}_{G,2,\varepsilon} \mathcal{T}_{B,2,\varepsilon} \hat{C}^{\dagger} \mathcal{T}_{G,1,\varepsilon} \mathcal{T}_{B,1,\varepsilon} \hat{C}^{\dagger} \mathcal{T}_{G,3,\varepsilon} \mathcal{T}_{B,3,\varepsilon}. \qquad (35) $$
It is easy to check that the corresponding continuum limit
is the same as in Eq. (34), without the factor 3 in front
of the first term on the rhs. The latter definition might be
more economical when it comes to implementing the QW on
experimental devices. However, we preferred to use the first
definition of $W_\epsilon$ because it can be more elegantly integrated
on the action of each of the $\tilde{T}_{B,\mu,\epsilon}$ operators, instead of an
isolated operation that only appears at the end of each time
step.

C. Spin network

In this work we have considered a particular simplicial
complex, namely a regular lattice of tetrahedra. It is quite
common to consider the dual graph of a simplicial complex,
namely by allocating, to each simplex a node, and to each
glueing between facets, a link between the nodes, see Fig. 7.
When the simplices are tetrahedra as in our case, the dual
graph is a four-valent graph quite close to a spin network, see
Fig. 8.

Spin networks form a basis of the Hilbert space of 3D
geometries in Loop Quantum Gravity. The questions of how
to propagate matter, such as massive Dirac particles, is a
long-standing problem in Loop Quantum Gravity. Our four-
valent graphs are variants, the main difference is that the
neighbours of each node are ordered, and that there are left-
handed and right-handed nodes, but those extra ingredients
seem crucial in order to achieve matter propagation. Indeed,
the present paper has just achieved a scheme which, upon the
spin network-like graph that is closest to flat space, reproduces
the Dirac equation in flat space. We wonder whether there
is a sense in which this exact same scheme, generalized to a
spin network (or a superposition of them) representing curved
space, reproduces the Dirac equation in this curved space.

<figure><img src="image_8.png" /><figcaption>FIG. 7: Simplicial complex and dual representation of a LH tetrahedron.</figcaption></figure>

V. SUMMARY AND PERSPECTIVES

In this paper we tessellate the 3D space with just two kinds
of tetrahedra, whose shapes are mirror of each other. We
place the amplitudes on their facets, and we alternate two
steps of dynamics according to Eq. (32): (i) we do a causal

shift, displacing amplitudes within the tetrahedron and across
their facets (ii) we apply a weakly localized mass coin, we
do a partial shift and apply a basis change coin, both strictly
localized to each tetrahedron. We prove that, in the continuum
limit, we obtain the Dirac equation in (3+1)-dimensions. By
taking the graph that is dual to the simplicial complex, we
obtain a variant of a Loop Quantum Gravity spin network.
The main difference is that the neighbors of each node
are ordered, and that there are left-handed and right-handed
nodes. Our scheme then provides a way to propagate spin
1/2 particles over these modified spin networks. Interestingly,
these modifications seem necessary in order to achieve such
matter propagation.

There are many possible avenues to explore. 3 + 1
Dirac QW on the grid have been extended to include the
electromagnetic field, or to account for the multi-particle
sector (quantum cellular automata). One could try to do
the same for tetrahedra. More interestingly perhaps, there
exists 3 + 1 ‘curved’ Dirac QW on the grid which, in the
continuum limit, yield the Dirac equation in curved spacetime,
as modelled by a metric field, which is used to parametrize
the coin. Alternatively, curved space can be modelled in
a combinatorial manner just by glueing tetrahedra in a way
that matches the discretization of 3D manifold. We wonder
whether running our Tetrahedral QW over some discretized
manifold, yields the same behaviour as the 3 + 1 curved Dirac
QW over the manifold, in some limit.

VI. ACKNOWLEDGEMENTS

This project/publication was made possible through the
support of the ID# 62312 grant from the John Templeton
Foundation, as part of the ‘The Quantum Information
Structure of Spacetime’ Project (QISS) . The opinions
expressed in this project/publication are those of the author(s)
and do not necessarily reflect the views of the John
Templeton Foundation. This work is supported by the
PEPR integrated project EPiQ ANR-22-PETQ-0007 part of
Plan France 2030, by the ANR JCJC DisQC ANR-22-
CE47-0002-01 founded from the French National Research
Agency and with the support of the french government under
the France 2030 investment plan, as part of the Initiative
d’Excellence d’Aix-Marseille Université - A*MIDEX AMX-
21-RID-011. From A.P. side, this work has been
founded by the Spanish MCIN/AEI/10.13039/501100011033
grant PID2020-113334GB-I00, Generalitat Valenciana grant
CIPROM/2022/66, the Ministry of Economic Affairs and
Digital Transformation of the Spanish Government through
the QUANTUM ENIA project call - QUANTUM SPAIN
project, and by the European Union through the Recovery,
Transformation and Resilience Plan - NextGenerationEU
within the framework of the Digital Spain 2026 Agenda,
and by the CSIC Interdisciplinary Thematic Platform (PTI+)
on Quantum Technologies (PTI-QTEP+). This project has
also received funding from the European Union’s Horizon
2020 research and innovation program under grant agreement
CaLIGOLA MSCA-2021-SE-01-101086123.
<figure><img src="image_9.png" /></figure>

(a) Causality of the black shift, represented with the arrows, which is the first step of the shifting dynamics. Arrows coming out of or going into the observed tetrahedra have a gradient effect to show the need to chained them up.

<figure><img src="image_10.png" /></figure>

(b) Strict locality of the grey shift and basis coin operator, the shift corresponds to the grey arrows (second step of the shifting dynamics) and the action of the strictly localized coin $\tilde{C}$ is represented with the blue and red lines. Weak locality of $\tilde{M}$ is shown by the dashed blue and red lines.

FIG. 8: Part of the spin network associated to the Dirac tetrahedral QW. Green and yellow nodes corresponds to the LH and RH tetrahedra, and two nodes are linked if their corresponding tetrahedra share a facet. The amplitudes lie on the edges. Therefore, a cube is represented by an hexagon. The position of the cubes in relation to each other is indicated by the label in their center. We show the action of the shift and coin operators for the two tetrahedra framed in black. In Fig. 8a we focus on the causal operations (the black shift), and in Fig. 8b we show the localized ones (grey shift, basis coin and mass coin).

## Appendix A: Same path

To show that the second walker follows the steps of the first one, we may compare the coordinates after one time step of a blue or red component to that of a cyan or magenta components initially located in the same tetrahedron. By looking at the direction of the arrows of Fig. 4b we see that in the LH (resp. RH) tetrahedra the direction of motion is positive (resp. negative) with respect to the three axis $x, y, z$ shown in Fig. 4a. Let $k_1(t)$ (resp. $k_2(t)$) be the coordinate of a blue or red (resp. cyan or magenta) component of the tetrahedron of coordinate $k$ at time $t$:

$$k_1(t) = \begin{pmatrix} k_x \\ k_y \\ k_z \end{pmatrix}, k_2(t) = k_1(t) \pm \varepsilon \begin{pmatrix} \kappa_x \\ \kappa_y \\ \kappa_z \end{pmatrix} \quad (A1)$$

where one of the $\kappa_\mu$ is equal to $1/2$ and the rest are 0. The non-vanishing element is located on the axis on which the first walk receives information at time $t$, as shown in Fig. 4a. After one time step, since a movement along one of the axis by a distance $\varepsilon$ has been made, it is straightforward noticing that:

$$k_1(t + \varepsilon) = k_1(t) \pm \varepsilon \begin{pmatrix} \kappa_x \\ \kappa_y \\ \kappa_z \end{pmatrix}, \quad (A2)$$

where one of the $\kappa_\mu$ is equal to 1 and the values of the others are identically 0. Similarly, by looking at Fig. 4b we see that:

$$k_2(t + \varepsilon) = k_2(t) \pm \varepsilon \begin{pmatrix} \kappa_x \\ \kappa_y \\ \kappa_z \end{pmatrix}. \quad (A3)$$

Now, two scenarios may arise. Either two of the $\kappa_\mu$ are $1/2$ and the remainder is 1, or the component is simply swaped on the other side of the facet on which it was, and so without involving a change of coordinates. There is only a difference of $\pm \varepsilon/2$ in one of the coordinates of $k_1$ and $k_2$ at time $t$. This proves that the second walker follows the first very closely. We show on Fig. 9 that the amplitudes of first and second walkers are following each others. Thus, we consider that they follow the same cubic path.

## Appendix B: Robust QW

As explained in the core of the paper, one of the main motivations of this study of QW over tetrahedra is to then be able to apply the same QW scheme over arbitrary graphs of degree bounded by four: whether to represent curvature in the combinatorial language of simplicial complexes; exotic crystallising structures; boundaries or defects. For instance in
<figure><img src="image_11.png" /><figcaption>FIG. 9: Path followed by amplitudes of the first and second walkers for the LH and RH tetrahedra after 3 time steps. The blue and red (resp. cyan and magenta) arrows are displayed with filled (resp. dashed) linestyle.</figcaption></figure>

[31, 32] defects are coded as missing nodes or broken links and located through a QW search.

Unfortunately the QW presented in the core of the paper is not robust to the presence of broken links, as it stands. To understand where the problem lies, comparison with the simpler 2 spatial dimensions scenario of Fig. 5 helps. Ignore the grey arrows and suppose that one of the nodes is missing. This leaves two black arrows dangling, with nowhere to go. The amplitudes that they carry cannot be dropped as this would break unitarity. Nor can they be placed back in their originating nodes, as new amplitudes have come to occupy the space.

The only possible fix is to double up the number of amplitudes, and balance the flow of information by considering the grey arrows. It seems unlikely that we can achieve robustness without doubling up the number of amplitudes. In the absence of a node, we shall be able to plug back the dangling outcoming black arrow of a node into its dangling incoming grey arrow. In this appendix, we do the Tetrahedral QW counterpart of this construction, in order to reach a less optimised, but more robust scheme.

# 1. Dynamics

As explained above, we can define a robust QW by doubling the number of amplitudes as shown in Fig. 10, and modifying the dynamics. We now only use swap operations for the shifting dynamics. This robust QW contains two independent walkers. The mirror walker does the same displacements as the first walker but in the opposite direction. Each facet of a tetrahedron owns one amplitude of each walker. Recall that the two-dimensional spinors lie in the facet shared by neighboring tetrahedra and contain an up and down component. Hence, each facet of a tetrahedron owns an up and down component of the two different walkers. The state of tetrahedron k at time t reads:

$$ \phi(t, k) = (\phi(t, k, i, \uparrow), \phi(t, k, i, \downarrow))_{i=0,...,3}^{T} \quad (B1) $$

The blue, red, cyan, magenta dots of Fig. 10 correspond to $i = 0, 1, 2, 3$. In the shifting dynamics of the non-robust first walker, the amplitudes located on the LH (resp. RH) tetrahedra move in a positive (resp. negative) direction over time. Moreover, we know from Eq. (7) that the up (resp. down) component of the non-robust first walker goes in negative (resp. positive) direction over time. Therefore, in the LH (resp. RH) tetrahedra, the up (resp. down) component is that of the mirror (resp. first) walker. Lastly, we define the j-th component of a tetrahedron as:

$$ \phi_j(t, k)_{j=0,...,7} = \begin{cases} \phi(t, k, j/2, \uparrow) & \text{if } j \text{ is even} \\ \phi(t, k, (j-1)/2, \downarrow) & \text{otherwise} \end{cases} \quad (B2) $$

## a. Shifting operator

The shifting dynamics is done in 3 steps. Inspecting Fig. 11, we see that for the LH (resp. RH) tetrahedra it consists in: (i) swap the mirror (resp. first) blue and the first (resp. mirror) cyan components, swap the first (resp. mirror) blue and mirror (resp. first) red components, swap the first (resp. mirror) red and mirror (resp. first) magenta components, (ii) swap the first and mirror blue components, swap the first and mirror red components, swap the first (resp. mirror) cyan and neighboring mirror (resp. first) magenta, swap the mirror (resp. first) magenta and the neighboring first (resp. mirror) cyan, (iii) swap the blue/cyan first (resp. mirror) components and swap the red/magenta mirror (resp. first) components. The time evolution of a tetrahedron reads:

$$ \phi(t + \varepsilon, k) = (S_2 S_1 S_0 \phi)(t, k) \quad (B3) $$

where the $S_i$ operators formally act on a tetrahedron as:
<figure><img src="image_12.png" /><figcaption>FIG. 10: Dissection of a cube into six tetrahedra for the robust QW. The big (resp. small) dots are that of the first (resp. mirror) walker.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>FIG. 11: Shifting dynamics of the robust tetrahedral QW. The black (filled), orange (thin) and pink (dotted) arrows correspond to the first, second and third step needed to implement our shifting dynamics. We have labelled the amplitudes for readability, after the several displacements the amplitudes moved from: 1 → 3, 2 → B, 4 → 2. Note that only the dynamics of the amplitudes of the middle tetrahedron is shown, thus, unrepresented displacements contain: A → 1, D → 4. Also, notice that only the coin step shuffles the amplitudes between LH and RH tetrahedra.</figcaption></figure>

$$ (S_0\phi)(t, k) = \begin{pmatrix} \phi_5(t, k) \\ \phi_2(t, k) \\ \phi_1(t, k) \\ \phi_6(t, k) \\ \phi_4(t, k) \\ \phi_0(t, k) \\ \phi_3(t, k) \\ \phi_7(t, k) \end{pmatrix}, \qquad (\text{B4a}) $$

$$ (S_1\phi)(t, k) = \begin{pmatrix} \phi_1(t, k) \\ \phi_0(t, k) \\ \phi_3(t, k) \\ \phi_2(t, k) \\ \phi_4(t, k) \\ \phi_6(t, n(k, 2)) \\ \phi_5(t, n(k, 3)) \\ \phi_7(t, k) \end{pmatrix}, \qquad (\text{B4b}) $$

$$ (S_2\phi)(t, k) = \begin{pmatrix} \phi_0(t, k) \\ \phi_5(t, k) \\ \phi_6(t, k) \\ \phi_3(t, k) \\ \phi_4(t, k) \\ \phi_1(t, k) \\ \phi_2(t, k) \\ \phi_7(t, k) \end{pmatrix}. \qquad (\text{B4c}) $$

Let us compute the state of a tetrahedron after one time step:
$$ (S_2 S_1 S_0 \phi)(t, k) = S_2 S_1 S_0 \begin{pmatrix} \phi_0(t, k) \\ \phi_1(t, k) \\ \phi_2(t, k) \\ \phi_3(t, k) \\ \phi_4(t, k) \\ \phi_5(t, k) \\ \phi_6(t, k) \\ \phi_7(t, k) \end{pmatrix} \\ = S_2 S_1 \begin{pmatrix} \phi_5(t, k) \\ \phi_2(t, k) \\ \phi_1(t, k) \\ \phi_6(t, k) \\ \phi_4(t, k) \\ \phi_0(t, k) \\ \phi_3(t, k) \\ \phi_7(t, k) \end{pmatrix} \tag{B5} \\ = S_2 \begin{pmatrix} \phi_2(t, k) \\ \phi_5(t, k) \\ \phi_6(t, k) \\ \phi_1(t, k) \\ \phi_4(t, k) \\ \phi_6(t, n(k, 2)) \\ \phi_5(t, n(k, 3)) \\ \phi_7(t, k) \end{pmatrix} \\ = \begin{pmatrix} \phi_2(t, k) \\ \phi_6(t, n(k, 2)) \\ \phi_5(t, n(k, 3)) \\ \phi_1(t, k) \\ \phi_4(t, k) \\ \phi_5(t, k) \\ \phi_6(t, k) \\ \phi_7(t, k) \end{pmatrix}. $$

We see that after one time step, only the blue ($j = 0, 1$) and red ($j = 2, 3$) components moved, while the others are left unchanged. Therefore, the global dynamics of our robust tetrahedral QW does not involve the cyan ($j = 4, 5$) and magenta ($j = 6, 7$) components. Thus, we only focus on the blue and red components. Although the cyan and magenta components do not move in the robust scheme after a time step, they are still necessary as ancillas in order to implement the shifting dynamics, as can be seen in Fig. 11. They act as intermediaries between adjacent tetrahedra to avoid instantaneous displacements. However, after one time step, these components have not actually moved. Therefore, they are simply a tool contributing to the smooth running of the robust displacement dynamics. In the non-robust scheme, the cyan and magenta components were moving, which is not the case in the robust scheme. Restricting shifting operations to swaps makes it difficult to replicate the non-robust dynamics where all components move. Moreover, the robust dynamics that comprise two walkers moving in opposite directions is more intuitive and straightforward than the walkers that follow each other of the non-robust scheme. Thus, it is more convenient to consider a simple dynamics involving two walkers moving in opposite directions, at the cost of

<figure><img src="image_14.png" /><figcaption>FIG. 12: The upper (resp. lower) figure show the global dynamics of the first (resp. mirror) walker. The global shifting dynamics of the robust tetrahedral QW is just the fusion of the two figures.</figcaption></figure>

immobilizing the cyan and magenta components, rather than trying to replicate the previous, less intuitive dynamic at all costs.

b. Coin operators

The coin operator of the first walker is $C = e^{i\frac{\pi}{3}} R_{\sigma_z}(\theta) R_{\sigma_x}(\theta)$ with $\theta = \pi/2$.

The mirror walker is doing the first walker walk in reverse. Thus, when applying the coin operator on the facets we want to map $\sigma_\mu$ into $\sigma_{\mu-1}$. Therefore, the coin of the mirror walker is $\tilde{C} = C^\dagger$.

In the following sections we will construct a robust QW that has the massive Dirac equation as its continuum limit. First, we will show how a Weyl equation can be simulated in the continuum limit. By combining two of such equations, we then can construct a robust QW for the massless Dirac equation. Finally, we will introduce the mass term to prevent travelling at the speed of light.
## 2. Weyl QW

### a. First walker

Each facet of a tetrahedron hosts a two dimensional spinor, let $\psi(t, k)$ (resp. $\tilde{\psi}(t, k)$) be that of the first (resp. mirror) walker. We formally define them as:

$$ \psi(t, k) = \begin{cases} \begin{pmatrix} \phi_0(t, n(k, 0)) \\ \phi_1(t, k) \\ \phi_2(t, k) \\ \phi_3(t, n(k, 1)) \end{pmatrix} & \text{if } k \text{ is LH} \\ \text{otherwise} & \end{cases}, \quad (B6) $$

$$ \tilde{\psi}(t, k) = \begin{cases} \begin{pmatrix} \phi_0(t, k) \\ \phi_1(t, n(k, 0)) \\ \phi_2(t, n(k, 1)) \\ \phi_3(t, k) \end{pmatrix} & \text{if } k \text{ is LH} \\ \text{otherwise} & \end{cases}. \quad (B7) $$

Let us first derive the continuum limit of the first walker. Interestingly, the first walker was constructed to have the Weyl equation as its continuum limit. Recall the Weyl equation can be rewritten as:

$$ i\partial_0\psi = D\psi \quad \text{with} \\ D = i\sum_j \sigma_j\partial_j. \quad (B8) $$

It is known that a translation by $\varepsilon$ along the $\mu$-axis approximates as:

$$ \tau_{\mu, \varepsilon} \psi = (\mathbb{I} + \varepsilon \partial_{\mu}) \psi + O(\varepsilon^2). \quad (B9) $$

Using the definition of $\sigma_3$ and Eq. (12) we have:

$$ (\mathbb{I} + \varepsilon\sigma_3\partial_3) \psi = T_{3,\varepsilon}\psi + O(\varepsilon^2). \quad (B10) $$

Given that $\sigma_1 = C\sigma_3 C^\dagger$, we also have:

$$ (\mathbb{I} + \varepsilon\sigma_1\partial_1) \psi = C T_{1,\varepsilon} C^\dagger \psi + O(\varepsilon^2). \quad (B11) $$

Finally along the $y$-axis, as $\sigma_2 = C\sigma_1 C^\dagger = CC\sigma_3 C^\dagger C^\dagger$ and $C^2 = C^\dagger$, we obtain:

$$ (\mathbb{I} + \varepsilon\sigma_2\partial_2) \psi = C^\dagger T_{2,\varepsilon} C\psi + O(\varepsilon^2). \quad (B12) $$

Lastly, the time evolution of the first walker is:

$$ \psi(t + \varepsilon, k) = (W_{\varepsilon}\psi)(t, k) + O(\varepsilon^2), \quad (B13) $$

with:

$$ \begin{aligned} W_{\varepsilon} &= C^{\dagger} T_{2,\varepsilon} C C T_{1,\varepsilon} C^{\dagger} T_{3,\varepsilon} \\ &= C^{\dagger} T_{2,\varepsilon} C^{\dagger} T_{1,\varepsilon} C^{\dagger} T_{3,\varepsilon}. \end{aligned} \quad (B14) $$

Using Eq. (B8) and Eq. (B9) we have:

$$ \psi(t + \varepsilon, k) = (\mathbb{I} + \varepsilon\partial_0)\psi(t, k) + O(\varepsilon^2). \quad (B15) $$

Moreover, with Eq. (B13) we get:

$$ \psi(t + \varepsilon, k) = \prod_{j=1}^{3} C^{\dagger}(\mathbb{I} + \varepsilon\sigma_3\partial_j)\psi(t, k) + O(\varepsilon^2) $$

$$ \psi(t + \varepsilon, k) = \prod_{j=1}^{3} (\mathbb{I} + \varepsilon\sigma_j\partial_j)\psi(t, k) + O(\varepsilon^2) $$

$$ (\mathbb{I} + \varepsilon\partial_0)\psi(t, k) = \left(\mathbb{I} + \varepsilon\sum_{j=1}^{3} \sigma_j\partial_j\right)\psi(t, k) + O(\varepsilon^2) $$

$$ \partial_0\psi(t, k) = \sum_{j=1}^{3} \sigma_j\partial_j\psi(t, k) + O(\varepsilon^2). \quad (B16) $$

In the continuum limit, we approximate the Weyl equation in (3+1)-dimensions, that is Eq. (B8).

### b. Mirror walker

We now focus on the mirror walker. The Weyl equation for the second chiral component has the form:

$$ i\partial_0\tilde{\psi} = \tilde{D}\tilde{\psi} \quad \text{with} \\ \tilde{D} = -i\sum_j \sigma_j\partial_j. \quad (B17) $$

Doing the same computations as before, we get for the $z$-axis:

$$ (\mathbb{I} - \varepsilon\sigma_3\partial_3)\tilde{\psi} = T_{3,\varepsilon}\tilde{\psi} + O(\varepsilon^2). \quad (B18) $$

Since $\sigma_1 = C\sigma_3 C^\dagger$, we also have:

$$ (\mathbb{I} - \varepsilon\sigma_1\partial_1)\tilde{\psi} = C T_{1,\varepsilon} C^\dagger \tilde{\psi} + O(\varepsilon^2). \quad (B19) $$

Lastly as $\sigma_2 = C^\dagger\sigma_3 C$, we obtain:

$$ (\mathbb{I} - \varepsilon\sigma_2\partial_2)\tilde{\psi} = C^\dagger T_{2,\varepsilon} C\tilde{\psi} + O(\varepsilon^2). \quad (B20) $$

The time evolution of the mirror walker is:

$$ \tilde{\psi}(t + \varepsilon, k) = \tilde{W}_{\varepsilon}\tilde{\psi}(t, k), \quad (B21) $$

with:

$$ \begin{aligned} \tilde{W}_{\varepsilon} &= C T_{1,\varepsilon} C^{\dagger} C^{\dagger} T_{2,\varepsilon} C T_{3,\varepsilon} \\ &= C T_{1,\varepsilon} C T_{2,\varepsilon} C T_{3,\varepsilon}. \end{aligned} \quad (B22) $$

In the continuum limit we do obtain the seeked Weyl
equation in (3+1)-dimensions since:

$$
\begin{align}
\tilde{\psi}(t + \varepsilon, k) &= \prod_{j=1}^{3} C(\mathbb{I} - \varepsilon\sigma_3 \partial_j)\tilde{\psi}(t, k) + O(\varepsilon^2) \nonumber \\
\tilde{\psi}(t + \varepsilon, k) &= \prod_{j=1}^{3} (\mathbb{I} - \varepsilon\sigma_j \partial_j)\tilde{\psi}(t, k) + O(\varepsilon^2) \nonumber \\
(\mathbb{I} + \varepsilon\partial_0)\tilde{\psi}(t, k) &= \left(\mathbb{I} - \varepsilon\sum_{j=1}^{3} \sigma_j \partial_j\right)\tilde{\psi}(t, k) + O(\varepsilon^2) \nonumber \\
\partial_0 \tilde{\psi}(t, k) &= -\sum_{j=1}^{3} \sigma_j \partial_j \tilde{\psi}(t, k) + O(\varepsilon^2). \tag{B23}
\end{align}
$$

3. Dirac QW

a. Massless Dirac QW

The four-dimensional spinor Ψ of the massless Dirac QW
contains all the amplitudes own by a facet where the first and
second (resp. third and fourth) components are that of the first
 resp. mirror) walker, that is Ψ(t, k) = ψ(t, k) ⊕ ˜ψ(t, k). In
(3+1)-dimensions, the massless Dirac equation is:

$$
\begin{align}
i \partial_0 \Psi &= \hat{D} \Psi \quad \text{with} \nonumber \\
\hat{D} &= \sigma_3 \otimes D \nonumber \\
&= i \sum_{j=1}^{3} \alpha^j \partial_j, \tag{B24}
\end{align}
$$

where D is that of Eq. (B8) and α^j = σ_3 ⊗ σ_j. The coin
operator is Ĉ = σ_0 ⊗ C, as α^{j+1} = Ĉα^j Ĉ†. Since the
dynamics of our robust QW takes place in 3 steps, let T_{μ,ε}^{(i)}
be the i-th step of translation of ε along the μ-axis. We start
with the z-axis:

$$
(\mathbb{I} + \varepsilon(\sigma_3 \otimes -\sigma_3)\partial_3) \Psi = \mathcal{T}_{3,\varepsilon}^{(2)} \mathcal{T}_{3,\varepsilon}^{(1)} \mathcal{T}_{3,\varepsilon}^{(0)} \Psi + O(\varepsilon^2). \quad (\text{B25})
$$

For the y-axis we get:

$$
(\mathbb{I} + \varepsilon(\sigma_3 \otimes -\sigma_2)\partial_2) \Psi = \hat{C}^\dagger \mathcal{T}_{2,\varepsilon}^{(2)} \mathcal{T}_{2,\varepsilon}^{(1)} \mathcal{T}_{2,\varepsilon}^{(0)} \hat{C} \Psi + O(\varepsilon^2). \quad (\text{B26})
$$

Lastly, for the x-axis we obtain:

$$
(\mathbb{I} + \varepsilon(\sigma_3 \otimes -\sigma_1)\partial_1) \Psi = \hat{C} \mathcal{T}_{1,\varepsilon}^{(2)} \mathcal{T}_{1,\varepsilon}^{(1)} \mathcal{T}_{1,\varepsilon}^{(0)} \hat{C}^\dagger \Psi + O(\varepsilon^2). \quad (\text{B27})
$$

The time evolution of the massless robust Dirac QW reads:

$$
\Psi(t + \varepsilon, k) = (\hat{W}_{\varepsilon}\Psi)(t, k), \quad (\text{B28})
$$

with:

$$
\hat{W}_{\varepsilon} = \prod_{j=1}^{3} \hat{C} \mathcal{T}_{j,\varepsilon}^{(2)} \mathcal{T}_{j,\varepsilon}^{(1)} \mathcal{T}_{j,\varepsilon}^{(0)}. \quad (\text{B29})
$$

We now verify that we approximate the Dirac equation:

$$
\begin{align*}
\Psi(t + \varepsilon, k) &= \prod_{j=1}^{3} \hat{C}(\mathbb{I} + \varepsilon(\sigma_3 \otimes -\sigma_3)\partial_j)\Psi(t, k) + O(\varepsilon^2) \\
\Psi(t + \varepsilon, k) &= \prod_{j=1}^{3} (\mathbb{I} + \varepsilon(\sigma_3 \otimes -\sigma_j)\partial_j)\Psi(t, k) + O(\varepsilon^2) \\
(\mathbb{I} + \varepsilon\partial_0)\Psi(t, k) &= \left(\mathbb{I} + \sum_{j=1}^{3} \varepsilon(\sigma_3 \otimes -\sigma_j)\partial_j\right)\Psi(t, k) + O(\varepsilon^2) \\
\partial_0\Psi(t, k) &= \sum_{j=1}^{3} (\sigma_3 \otimes -\sigma_j)\partial_j\Psi(t, k) + O(\varepsilon^2). \tag{B30}
\end{align*}
$$

In the continuum limit we approximate Eq. (B24).

b. Massive Dirac QW

We add the mass term to the Dirac QW in the same way we
did in Sec. IV B, namely by adding an additional mass coin.
Hence, the Dirac equation in Planck units reads:

$$
i \partial_0 \Psi = \mathcal{D} \Psi \quad \text{with} \\
\mathcal{D} = m \gamma^0 + i \sum_j \alpha^j \partial_j, \qquad (\text{B31})
$$

where *m* is the mass. The matrix γ⁰ is that of Eq. (29), i.e.

$$
\gamma^0 = \sigma_1 \otimes \sigma_0. \tag{B32}
$$

Let $M = e^{-i\epsilon m \gamma^0}$ be the associated mass coin. Moreover, we notice that

$$
\begin{align}
\hat{C} M \hat{C}^{\dagger} &= M \tag{B33} \\
M \hat{C} M^{\dagger} &= \hat{C}, \notag
\end{align}
$$

making it possible to apply both coins after translating along
each axis. Hence, the time evolution of the robust Dirac
tetrahedral QW is:

$$
\Psi(t + \varepsilon, k) = (W_{\varepsilon}\Psi)(t, k), \quad (\text{B34})
$$

where

$$
\hat{W}_{\varepsilon} = \prod_{j=1}^{3} \hat{C} \mathcal{T}_{j,\varepsilon}^{(2)} \mathcal{T}_{j,\varepsilon}^{(1)} \mathcal{T}_{j,\varepsilon}^{(0)} M. \quad (\text{B35})
$$

Lastly, in the continuum limit we get:

$$
\begin{align*}
\Psi(t + \varepsilon, k) &= (\mathbb{I} - i\varepsilon\mathcal{D})\Psi(t, k) + O(\varepsilon^2) \\
&= \prod_{j=1}^{3} \hat{C}^{\dagger}(\mathbb{I} + \varepsilon(\sigma_3 \otimes \sigma_3)\partial_j)M\Psi(t, k) + O(\varepsilon^2) \\
&= \prod_{j=1}^{3} (\mathbb{I} + \varepsilon(\sigma_3 \otimes \sigma_j)\partial_j)M\Psi(t, k) + O(\varepsilon^2).
\tag{B36}
\end{align*}
$$
Moreover, as $M = \mathbb{I} - i\epsilon m\gamma^0 + O(\epsilon^2)$, we have:

$$ (\mathbb{I} + \epsilon \partial_0) \Psi(t, k) = \left( \mathbb{I} - 3im\epsilon\gamma^0 + \sum_{j=1}^{3} \epsilon(\sigma_3 \otimes \sigma_j)\partial_j \right) \Psi(t, k) + O(\epsilon^2) $$

$$ \partial_0 \Psi(t, k) = \left( -3im\gamma^0 + \sum_{j=1}^{3} (\sigma_3 \otimes \sigma_j) \partial_j \right) \Psi(t, k) + O(\epsilon^2), \quad (B37) $$

[1] Andris Ambainis. Quantum walk algorithm for element distinctness. *SIAM Journal on Computing*, 37(1):210–239, 2007. I

[2] Quentin Aristote, Nathanaël Eon, and Giuseppe Di Molfetta. Dynamical triangulation induced by quantum walk. *Symmetry*, 12(1):128, 2020. I

[3] Pablo Arnault and Fabrice Debbasch. Quantum walks and discrete gauge theories. *Physical Review A*, 93(5):052301, 2016. I

[4] Pablo Arnault and Fabrice Debbasch. Quantum walks and gravitational waves. *Annals of Physics*, 383:645–661, 2017. I

[5] Pablo Arnault, Benjamin Pepper, and A Pérez. Quantum walks in weak electric fields and bloch oscillations. *Physical Review A*, 101(6):062324, 2020. I

[6] Pablo Arrighi, Giuseppe Di Molfetta, Iván Márquez-Martín, and Armando Pérez. Dirac equation as a quantum walk over the honeycomb and triangular lattices. *Physical Review A*, 97(6):062111, 2018. I

[7] Pablo Arrighi, Giuseppe Di Molfetta, Ivan Marquez-Martín, and Armando Perez. From curved spacetime to spacetime-dependent local unitaries over the honeycomb and triangular quantum walks. *Scientific Reports*, 9(1):10904, July 2019. I

[8] Pablo Arrighi and Stefano Facchini. Quantum walking in curved spacetime: (3 + 1) dimensions, and beyond, 2016. I

[9] Pablo Arrighi, Stefano Facchini, and Marcelo Forets. Discrete lorentz covariance for quantum walks and quantum cellular automata. *New Journal of Physics*, 16(9):093007, 2014. I

[10] Pablo Arrighi, Vincent Nesme, and Marcelo Forets. The dirac equation as a quantum walk: higher dimensions, observational convergence. *Journal of Physics A: Mathematical and Theoretical*, 47(46):465302, 2014. I

[11] János K Asbóth. Symmetries, topological phases, and bound states in the one-dimensional quantum walk. *Physical Review B*, 86(19):195414, 2012. I

[12] Andrea Barbieri. Quantum tetrahedra and simplicial spin networks. *Nuclear Physics B*, 518(3):714–728, 1998. I

[13] Khadeejah Bepari, Sarah Malik, Michael Spannowsky, and Simon Williams. Quantum walk approach to simulating parton showers. *Physical Review D*, 106(5):056002, 2022. I

[14] Eugenio Bianchi, Muxin Han, Carlo Rovelli, Wolfgang Wieland, Elena Magliaro, and Claudio Perini. Spininfoam fermions. *Classical and Quantum Gravity*, 30(23):235023, oct 2013. I

[15] CM Chandrashekar. Two-component dirac-like hamiltonian for

which approximates Eq. (B31).

Finally, the robust Dirac tetrahedral QW is summarized in Eq. (B35). It consists of: (i) applying a mass coin in a strictly localized manner, (ii) making three translations, where the second one is strictly localized, and the first and last ones are weakly localized, (iii) applying a directional coin that is weakly localized to the tetrahedra. The different stages can not be merged due to their slightly different locality. Indeed, even though the directional coin $\hat{C}$ and the second translation step $\mathcal{T}_{j,\epsilon}^{(1)}$ are strictly local, and the remaining translations $\mathcal{T}_{j,\epsilon}^{(0,2)}$ with the mass coin $M$ are weakly local, the merger is impossible, as the directional coin must be applied after the translations, and these do not commute.

generating quantum walk on one-, two- and three-dimensional lattices. *Scientific reports*, 3(1):2829, 2013. I

[16] Francesco De Nicola, Linda Sansoni, Andrea Crespi, Roberta Ramponi, Roberto Osellame, Vittorio Giovannetti, Rosario Fazio, Paolo Mataloni, and Fabio Sciarrino. Quantum simulation of bosonic-fermionic noninteracting particles in disordered systems via a quantum walk. *Physical Review A*, 89(3):032322, 2014. I

[17] Jean-Yves Desaules, Kieran Bull, Aiden Daniel, and Zlatko Papić. Hypergrid subgraphs and the origin of scarred quantum walks in many-body hilbert space. *Physical Review B*, 105(24):245137, 2022. I

[18] Giuseppe Di Molfetta, Marc Brachet, and Fabrice Debbasch. Quantum walks as massless dirac fermions in curved spacetime. *Physical Review A*, 88(4):042301, 2013. I

[19] Giuseppe Di Molfetta, Marc Brachet, and Fabrice Debbasch. Quantum walks in artificial electric and gravitational fields. *Physica A: Statistical Mechanics and its Applications*, 397:157–168, 2014. I

[20] Giuseppe Di Molfetta and Fabrice Debbasch. Discrete-time quantum walks: Continuous limit and symmetries. *Journal of Mathematical Physics*, 53(12), 2012. I

[21] Nathanaël Eon, Giuseppe Di Molfetta, Giuseppe Magnifico, and Pablo Arrighi. A relativistic discrete spacetime formulation of 3+1 qed. *arXiv preprint arXiv:2205.03148*, 2022. I

[22] John King Gamble, Mark Friesen, Dong Zhou, Robert Joynt, and SN Coppersmith. Two-particle quantum walks applied to the graph isomorphism problem. *Physical Review A*, 81(5):052313, 2010. I

[23] Salvador Godoy, Lucio Andrade, and E Braun. Quantum random walk and transient phenomena. *Physica A: Statistical Mechanics and its Applications*, 196(3):416–426, 1993. I

[24] Salvador Godoy and Shigeji Fujita. A quantum random-walk model for tunneling diffusion in a 1d lattice. a quantum correction to fick's law. *The Journal of Chemical Physics*, 97(7):5148–5154, 1992. I

[25] Gareth Jay, Fabrice Debbasch, and Jingbo B Wang. Dirac quantum walks on triangular and honeycomb lattices. *Physical Review A*, 99(3):032113, 2019. I

[26] Nicolas Jolly and Giuseppe Di Molfetta. Twisted quantum walks, generalised dirac equation and fermion doubling. *The European Physical Journal D*, 77(5):80, 2023. I

[27] Ivan Márquez-Martín, Giuseppe Di Molfetta, and Armando Pérez. Fermion confinement via quantum walks in (2+ 1)-dimensional and (3+ 1)-dimensional space-time. *Physical
Review A, 95(4):042112, 2017. I

[28] Palash B. Pal. Dirac, majorana, and weyl fermions. Am. J. Phys., 79(5):485-498, May 2011. III

[29] Michael E. Peskin and Daniel V. Schroeder. *An Introduction to quantum field theory*. Addison-Wesley, Reading, USA, 1995. III, III

[30] Vinay V Ramasesh, Emmanuel Flurin, M Rudner, Irfan Siddiqi, and Norman Y Yao. Direct probe of topological invariants using bloch oscillating quantum walks. *Physical review letters*, 118(13):130501, 2017. I

[31] Mathieu Roget, Stéphane Guillet, Pablo Arrighi, and Giuseppe Di Molfetta. Grover search as a naturally occurring phenomenon. *Physical Review Letters*, 124(18):180501, 2020. I, B

[32] Mathieu Roget, Hachem Kadri, and Giuseppe Di Molfetta. Optimality conditions for spatial search with multiple marked vertices. Phys. Rev. Res., 5:033021, Jul 2023. I, B

[33] Carlo Rovelli and Lee Smolin. Spin networks and quantum gravity. Physical Review D, 52(10):5743, 1995. I

[34] Muhammad Sajid, János K Asbóth, Dieter Meschede, Reinhard F Werner, and Andrea Alberti. Creating anomalous floquet chern insulators with magnetic quantum walks. Physical Review B, 99(21):214303, 2019. I

[35] Kevissen Sellapillay, Pablo Arrighi, and Giuseppe Di Molfetta. A discrete relativistic spacetime formalism for 1+ 1-qed with continuum limits. Scientific Reports, 12(1):2198, 2022. I

[36] Kevissen Sellapillay, Alberto D Verga, and Giuseppe Di Molfetta. Entanglement dynamics and ergodicity breaking in a quantum cellular automaton. Physical Review B, 106(10):104309, 2022. I

[37] Nicholas Slate, Edric Matwiejew, Samuel Marsh, and Jingbo B Wang. Quantum walk-based portfolio optimisation. Quantum, 5:513, 2021. I

[38] D. M. Y. Sommerville. Space-filling tetrahedra in euclidean space. Proceedings of the Edinburgh Mathematical Society, 41:49-57, 1922. I, II

[39] Frederick W Strauch. Relativistic quantum walks. Physical Review A, 73(5):054302, 2006. I

[40] Alberto Verga and Ricardo Gabriel Elias. Entanglement and interaction in a topological quantum walk. SciPost Physics, 5(2):019, 2018. I

[41] Julien Zylberman, Giuseppe Di Molfetta, Marc Brachet, Nuno F Loureiro, and Fabrice Debbasch. Quantum simulations of hydrodynamics via the madelung transformation. Physical Review A, 106(3):032408, 2022. I