Introduction

Quantum walks (QWs) are the quantum analog of classical random walks, in which the walker steps forward or backward along a line based on a coin flip. In a QW, the walker proceeds in a quantum superposition of paths, and the resulting interference forms the basis of a wide variety of quantum algorithms, such as quantum search1,2,3,4,5, graph isomorphism problems6,7,8, ranking nodes in a network9,10,11,12, and quantum simulations, which mimic different quantum systems at the low and high energy scale13,14,15,16,17,18,19,20,21,22. In the discrete-time QW (DQW)23,24, a quantum coin operation is introduced to prescribe the direction in which the particle moves in position space at each discrete step. In the continuous-time QW (CQW)25,26, one can directly define the walk evolution on position space itself using continuous-time evolution. We focus on DQWs and their implementation on gate-based quantum circuits in this work.

DQWs can be realized directly on lattice-based quantum systems where position space matches the discrete lattice sites. Such implementations have been reported with cold atoms27,28 and photonic systems29,30,31,32. In trapped ions, a DQW has been implemented by mapping position space to locations in phase space given by the degrees of freedom associated with the harmonic motion of the ion in the trap33,34,35. All these physical implementations have followed an analogue quantum simulation approach. However, implementing QWs on a circuit-based system is crucial to explore the algorithm applications based on QWs. The implementation of a DQW on a three-qubit NMR system36, a CQW on a two-qubit photonic processor37 and a split-step QW on superconducting circuits38,39 are the circuit-based implementations reported to date. To implement DQWs on circuit-based quantum processors, its necessary to map the position space to the available multi-qubit states. The range of the walk is set by the available qubit number and gate depth. The term Quantum Cellular Automaton (QCA) describes a unitary evolution of a particle on a discretized space40,41,42, as occurs with QWs. In this context, the one-dimensional Dirac cellular automaton (DCA) has been derived from the symmetries of the QCA showing how the dynamics of the Dirac equation emerges40,41,42,43,44.

Here we implement efficient quantum circuits for a DQW in one-dimensional position space, which provide the time-evolution up to five steps. We report the experimental realization of a DQW on five qubits within a seven-qubit programmable trapped-ion quantum computer45. With a tunable walk probability at each step we also show the experimental realization of a DCA where the coin bias parameter mimics the mass term in the Dirac equation. This will be central for discrete-time quantum simulation of the dynamics associated with the relativistic motion of a spin-1/2 particle in position space.

Results

Review of quantum walks and the connection to the Dirac equation

The DQW consists of two quantum mechanical systems, an effective coin and the position space of the walker, as well as an evolution operator, which is applied to both systems in discrete time-steps. The evolution is given by a unitary operator defined on a tensor product of two Hilbert spaces \({{\mathcal{H}}}_{{\rm{c}}}\otimes {{\mathcal{H}}}_{{\rm{p}}}\) where, \({{\mathcal{H}}}_{{\rm{c}}}\) is the coin Hilbert space spanned by the internal states \({\left|0\right\rangle }_{{\rm{c}}}\) and \({\left|1\right\rangle }_{{\rm{c}}}\) of a single qubit, while \({{\mathcal{H}}}_{p}\) represents the position Hilbert space given by the position states \(\left|x\right\rangle\) with \(x\in {\mathbb{Z}}\) encoded in several qubits as described below. Here, the unitary quantum coin toss operation, \({\hat{C}}_{\theta }\), is a unitary rotation operator that acts on the coin qubit space,

$${\hat{C}}_{\theta }=\left[\begin{array}{cc}\cos \theta &-i\sin \theta \\ -i\sin \theta &\cos \theta \end{array}\right]\otimes {\hat{I}}_{{\rm{p}}},$$
(1)

where θ is a coin bias parameter that can be varied at each step to modify the QW path superposition weights. The conditional position-shift operator, \(\hat{S}\), translates the particle to the left and right conditioned by the state of the coin qubit,

$$\hat{S}={\left|0\right\rangle }_{{\rm{c}}\,{\rm{c}}}\langle 0| \otimes \sum _{x\in {\mathbb{Z}}}| x-1\rangle {\langle x| +| 1\rangle }_{{\rm{c}}\,{\rm{c}}}\langle 1| \otimes \sum _{x\in {\mathbb{Z}}}| x+1\rangle \left\langle x\right|.$$
(2)

The state of the particle in position space after t steps of the walk, is accomplished by the repeated action of the operator \(\hat{W}=\hat{S}{\hat{C}}_{\theta }\) on the initial state of the particle \({\left|\psi \right\rangle }_{{\rm{c}}}=\alpha {\left|0\right\rangle }_{{\rm{c}}}+\beta {\left|1\right\rangle }_{{\rm{c}}}\) at position x = 0, as shown in Fig. 1,

$$\left|\Psi (x,t)\right\rangle ={\hat{W}}^{t}\left[{\left|\psi \right\rangle }_{c}\otimes \left|x=0\right\rangle \right]=\sum _{x}\left[\begin{array}{c}{\psi }_{x,t}^{0}\\ {\psi }_{x,t}^{1}\end{array}\right],$$
(3)

where \({\psi }_{x,t}^{0(1)}\) denotes the left(right) propagating component of the particle at time-step t. The probability of finding the particle at position x and time t will be \(P(x,t)=| {\psi }_{x,t}^{0}{| }^{2}+| {\psi }_{x,t}^{1}{| }^{2}\).

Fig. 1: Discrete-time quantum walk scheme.
figure 1

Each step is composed of a quantum coin operation, \({\hat{C}}_{\theta }\), with tunable effective coin bias parameters, θi, followed by a shift operation, \(\hat{S}\).

Recent works have shown a relationship between DQWs and the Dirac equation14,15,16,17,18,43. Starting form a discrete-time evolution operator and then moving from position space to momentum space, Dirac kinematics can be recovered from the diagonal terms of the unitary evolution operator for small momenta in the small mass regime16,17,18. In contrast with these proposals in the Fourier frame, we focus our implementation on the probability distribution of the DQW, which is analogous to the spreading of a relativistic particle. To realize a DCA and recover the Dirac equation, a split-step quantum walk, one form of the DQW, is used40. Each step of a split-step quantum walk is a composition of two half step evolutions with different coin biases and position-shift operators,

$${\hat{W}}_{{\rm{ss}}}={\hat{S}}_{+}{\hat{C}}_{{\theta }_{2}}{\hat{S}}_{-}{\hat{C}}_{{\theta }_{1}},$$
(4)

where the coin operation \({\hat{C}}_{{\theta }_{j}}\), with j = 1, 2, is given in Eq. (1). The split-step position-shift operators are,

$${\hat{S}}_{-}={\left|0\right\rangle }_{{\rm{c}}\,{\rm{c}}}\langle 0| \otimes \sum _{x\in {\mathbb{Z}}}| x-1\rangle {\langle x| +| 1\rangle }_{{\rm{c}}\,{\rm{c}}}\langle 1| \otimes \sum_{x\in {\mathbb{Z}}}| x\rangle \left\langle x\right|,$$
(5)
$${\hat{S}}_{+}={\left|0\right\rangle }_{{\rm{c}}\,{\rm{c}}}\langle 0| \otimes \sum_{x\in {\mathbb{Z}}}| x\rangle {\langle x| +| 1\rangle }_{{\rm{c}}\,{\rm{c}}}\langle 1| \otimes \sum_{x\in {\mathbb{Z}}}| x+1\rangle \left\langle x\right|.$$
(6)

Following Mallick40 and Kumar44, the particle state at time t and position x after the evolution operation \({\hat{W}}_{{\rm{ss}}}\) is described by the differential equation,

$$\frac{\partial }{\partial t}\left[\begin{array}{c}{\psi }_{x,t}^{0}\\ {\psi }_{x,t}^{1}\end{array}\right]= \cos {\theta }_{2}\left[\begin{array}{cc}\cos {\theta }_{1} & -i\sin {\theta }_{1}\\ i\sin {\theta }_{1} & -\cos {\theta }_{1}\end{array}\right]\left[\begin{array}{c}\frac{\partial {\psi }_{x,t}^{0}}{\partial x}\\ \frac{\partial {\psi }_{x,t}^{1}}{\partial x}\end{array}\right]\\ +\left[\begin{array}{cc}\cos ({\theta }_{1}+{\theta }_{2})-1 & -i\sin ({\theta }_{1}+{\theta }_{2})\\ -i\sin ({\theta }_{1}+{\theta }_{2}) & \cos ({\theta }_{1}+{\theta }_{2})-1\end{array}\right]\left[\begin{array}{c}{\psi }_{x,t}^{0}\\ {\psi }_{x,t}^{1}\end{array}\right].$$
(7)

The tunability of parameters θ1 and θ2 on the split-step QW permits the study of one-dimensional Dirac equations effectively, within the low momentum subspace, for spin-1/2 particles40,44. It is important to stress out that, the description of the Dirac equation used here corresponds to the 2 × 2 representation, i.e. no spin degree of freedom. For instance, the massless particle Dirac equation can be recovered for \(\cos ({\theta }_{1}+{\theta }_{2})=1\). Thereby, Eq. (7) becomes \(i\hslash [{\partial }_{t}-\cos {\theta }_{2}(\cos {\theta }_{1}{\sigma }_{z}+\sin {\theta }_{1}{\sigma }_{y}){\partial }_{x}]\Psi (x,t)=0\), which is identical to the Dirac equation of a massless particle in the relativistic limit46. In contrast, considering θ1 = 0 and a very small value of θ2 corresponds to the Dirac equation for particles with small mass35,46 in the form \(i\hslash [{\partial }_{t}-(1-{\theta }_{2}^{2}/2){\sigma }_{z}{\partial }_{x}+i{\theta }_{2}{\sigma }_{x}]\Psi (x,t)\approx 0\).

At the same time, by choosing θ1 = 0, the quantum walk operator \({\hat{W}}_{{\rm{ss}}}\) given in Eq. (4) takes the form of the unitary operator for a DCA40,

$${\hat{W}}_{{\rm{ss}}}=\left[\begin{array}{cc}\cos ({\theta }_{2}){S}_{-}&-i\sin ({\theta }_{2}){\mathbb{1}}\\ -i\sin ({\theta }_{2}){\mathbb{1}}&\cos ({\theta }_{2}){S}_{+}\end{array}\right]={U}_{{\rm{DCA}}}.$$
(8)

Within this framework, θ2 determines the mass of the Dirac particle. The split-step DQW described by the operator \({\hat{W}}_{{\rm{ss}}}\) is equivalent to the two period DQW with alternate coin operations, θ1 and θ2, when the alternate points in position space with zero probability are ignored47. Therefore, all the dynamics of a DCA can be recovered from the DQW evolution using \(\hat{W}\) and alternating the two coin operations. See Methods for a comparison between DCA and the explicit solution of the Dirac equation. Typical features of the Dirac equation in relativistic quantum mechanics, such as the Zitterbewegung40 and the Klein paradox48, are also dynamical features of the DCA, as well as the spreading of the probability distribution and the entanglement of localized positive-energy states. We note that these effects have also been shown in direct analog simulations of the Dirac equation with trapped ions35 and BECs49.

Experimental DQW implementation

To realize the DQW on a system of qubits one must pick a mapping of the particle position to the qubit space. As shown in50, there is no unique way to map position states to multi-qubit states, so each circuit decomposition depends on the configuration adopted. A direct mapping of each walker position to one qubit in the chain mimicking the arrangement of the qubit array is inefficient in terms of qubit number and gates required (the former grows linearly and the latter quadratically with the position space size modeled). In order to minimize resource use, we take advantage of a digital representation to map the position space into a multi-qubit state and re-order it in such a way that the state \(\left|0\right\rangle \,(\left|1\right\rangle )\) of the last qubit corresponds to even (odd) position numbers. This allows us to minimize the changes needed in the qubit space configuration during each step of the walk (see Fig. 2). To implement a quantum walk in one-dimensional position Hilbert space of size 2n, (n + 1) qubits are required. One qubit acts as the coin and the other n qubits mimic the position Hilbert space with 2n − 1 positions of a symmetric walk about \(\left|x=0\right\rangle\). We note that the particle can be started from any point in the position space, however setting the initial state reduces the gate counting in the circuit and hence reduces the overall error. The coin operation is achieved by single-qubit rotations on the coin-qubit while the shift operators are realized by using the coin as a control qubit to change the position state during the walk.

Fig. 2: Mapping of multi-qubit states to position states.
figure 2

Multi-qubit states are re-ordered in such a way that the state \(\left|0\right\rangle \,(\left|1\right\rangle )\) of the last qubit corresponds to even (odd) position numbers and its correspondence in the position space.

We realize the walk on a chain of seven individual 171Yb+ ions confined in a Paul trap and laser-cooled close to their motional ground state45,51. Five of these are used to encode qubits in their hyperfine-split 2S1/2 ground level. Single-qubit rotations, or R gates, and two-qubit entangling interactions, or XX gates are achieved by applying two counter-propagating optical Raman beams to the chain, one of which features individual addressing (see Methods for experimental details). We can represent up to 15 positions of a symmetric QW, including the initial position \(\left|x=0\right\rangle\).

Based on this position representation a circuit diagram for the DQW on five qubits with the initial state \({\left|0\right\rangle }_{c}\otimes \left|0000\right\rangle\) is composed for up to five steps, see Fig. 3. Each evolution step, \(\hat{W}\), starts with a rotation operation on the coin-qubit, \({\hat{C}}_{{\theta }_{j}}\), followed by a set of controlled gates that change the position state of the particle under \(\hat{S}\). Due to the gratuitous choice of position representation used, it is enough to perform a single-qubit rotation on the last qubit at every step, which could also be done by classical tracking50.

Fig. 3: Circuit implementation of quantum walks on a trapped-ion processor and its time evolution.
figure 3

a Circuit diagram for a DQW and DCA. Each dashed block describes one step in the quantum walk. b Discrete-time Quantum Walk. Comparison of the experimental results (left) and the theoretical quantum-walk probability distribution (right) for the first five steps with initial particle state b i and b iv \({\left|\psi \right\rangle }_{{\rm{c}}}={\left|0\right\rangle }_{{\rm{c}}}\), b ii and b iv \({\left|\psi \right\rangle }_{{\rm{c}}}={\left|1\right\rangle }_{{\rm{c}}}\), b iii and b vi \({\left|\psi \right\rangle }_{{\rm{c}}}={\left|0\right\rangle }_{{\rm{c}}}+i{\left|1\right\rangle }_{{\rm{c}}}\), and position state \(\left|x=0\right\rangle\). c Output of a step-5 Dirac Cellular Automaton for θ1 = 0 and, c i and c iv θ2 = π/4, c ii and c v θ2 = π/10 and c iii and c vi θ2 = π/20 with the initial state \(\left|{\Psi }_{{\rm{in}}}\right\rangle =({\left|0\right\rangle }_{{\rm{c}}}+i{\left|1\right\rangle }_{{\rm{c}}})\otimes \left|x=0\right\rangle\).

Computational gates such as CNOT, Toffoli, and Toffoli-4 are generated by a compiler which breaks them down into constituent physical-level single- and two-qubit gates45. A circuit diagram detailing the compiled building blocks is shown in Methods. To prepare an initial particle state different from \({\left|0\right\rangle }_{{\rm{c}}}\) it is enough to perform a rotation on the coin-qubit before the first step. In some cases this rotation can be absorbed into the first gates in step one. Table 1 summarizes the number of native gates needed per step for initial state. To recover the evolution of the Dirac equation in a DQW after five steps, 81 single qubit gates and 32 XX-gates are required.

Table 1 Gate counting.

After evolving a number of steps, we sample the corresponding probability distribution 3000 times and correct the results for readout errors. For the DQW evolution up to five steps shown in Fig. 3, a balanced coin (θ1 = θ2 = π/4) is used where the initial position is \(\left|x=0\right\rangle\) for different initial particle states, \({\left|0\right\rangle }_{{\rm{c}}}\) in Fig. 3b i, \({\left|1\right\rangle }_{{\rm{c}}}\) in Fig. 3b ii, and an equal superposition of both in Fig. 3b iii. In Fig. 3b iv, b v, and b vi we show the ideal output from classical simulation of the circuit for comparison (see Methods for a plot of the difference). With a balanced coin the particle evolves in equal superposition to the left and right position at each time step and upon measurement, there is a 50/50 probability of finding the particle to the left or right of its previous position, just as in classical walk. If we let the DQW evolve for more than three steps before we perform a position measurement, we will find a very different probability distribution compared to the classical random walk52.

The same experimental setup can be used to recover a DCA with a two-period DQW. Here we set θ1 = 0 and varied θ2 to recover the Dirac equation for different mass values. In Fig. 3c, we show experimental results for θ2 = π/4, π/10 and π/20, corresponding to a mass 1.1357, 0.3305, and 0.159 in units of c−2s−1, with the initial particle state in the superposition \({\left|0\right\rangle }_{c}+i{\left|1\right\rangle }_{c}\). The main signature of a DCA for small mass values is the presence of peaks moving outward and a flat distribution in the middle as shown for the cases with small values of θ2, Figs. 3c ii-iii. This bimodal probability distribution in position space is an indication of the one-dimensional analog of an initially localized Dirac particle, with positive energy, evolving in time which spreads in all directions in position space at speeds close to the speed of light53. In contrast, a DCA with θ2 = π/4, Fig. 3c i corresponds to a massive particle and hence there is a slow spread rather than a ballistic trajectory in position space.

Discussion

We have shown how quantum walks form the basic elements for simulation of the dynamics associated with the free Dirac particle with positive energy. Despite the population mismatch of 0.05–0.2 between the simulation and the experimental results after five steps, the final probability density exhibits the characteristic behavior of an initially localized Dirac particle. A key factor on the digitization of DQW/DCA is the mapping of qubit states to position space. An adequate mapping is important to minimize the number of gates on the protocol, and as a consequence, the resource scaling of the evolution. By increasing in the available number of qubits, these quantum circuits can be scaled to implement more steps and simulate a multi-particle DQW. The number of gates has a polynomial growth rate with the number of steps54. The correspondence between DQWs and the dynamics of Dirac particles suggests that the QWs formalism is as a viable approach to reproduce a variety of phenomena underpinned by Dirac particle dynamics in both the high- and low-energy regime22,39,43. Quantum simulations of free quantum field theory43, Yang-Mills gauge-field on fermionic matter55, as well as the effect of mass and space-time curvature on entanglement between accelerated particles20,56,57 have been reported and probing quantum field theory from the algorithmic perspective in an active field of research. However, the circuit complexity for position-dependent coin operations needed for simulating these effects will increase with the complexity of the evolution, which means further improvements in quantum hardware will be necessary for their realization.

Methods

Experimental details

The experiments are performed in a chain of seven individual 171Yb+ ions confined in a Paul trap and laser-cooled close to their motional ground state45,51. In order to guarantee higher uniformity in the ion spacing, matching the equally spaced individual addressing beams, the middle five of these are used to encode qubits in their hyperfine-split 2S1/2 ground level, with an energy difference of 12.642821 GHz. The two edge ions are neither manipulated nor measured, however, their contribution to the collective motion is included when creating the entangling operations. The ions are initialized by an optical pumping scheme and are collectively read out using state-dependent fluorescence detection58, with each ion being mapped to a distinct photo-multiplier tube (PMT) channel. The system has two mechanisms for quantum control, which can be combined to implement any desired operation: single-qubit rotations, or R gates, and two-qubit entangling interactions, or XX gates. These quantum operations are achieved by applying two counter-propagating optical Raman beams from a single 355-nm mode-locked laser59. The first Raman beam is a global beam applied to the entire chain, while the second is split into individual addressing beams, each of which can be controlled independently and targets one qubit. Single-qubit gates are generated by driving resonant Rabi rotations of defined phase, amplitude, and duration. Two-qubit gates are realized by illuminating two ions with beat-note frequencies near to the motional sidebands and creating an effective spin-spin (Ising) interaction via transient entanglement between the state of two ions and all modes of motion60,61,62. The average state detection fidelity for single- and two-qubit gate are 99.5(2)% and 98–99%, respectively. Rotations around the z-axis are achieved by phase advances on the classical control signals. Both the R as well as the XX angle can be varied continuously. State preparation and measurement (SPAM) errors are characterized and corrected by applying the inverse of an independently measured state-to-state error matrix63.

Errors

In order to illustrate how our experiment performs, we plot the absolute value of the difference between measured and simulated position distributions, Fig. 4, they match the theoretical expectation closely. These distributions are obtained after tracing out the coin information of the unitary evolution \(\hat{W}=\hat{S}{\hat{C}}_{\theta }\) for each time-step. In both instances, DQW and DCA, the number of gates and hence the error incurred grows with the number of steps.

Fig. 4: Experimental errors.
figure 4

Experimental error distribution for DQW (left) and DCA (right).

Apart from this, the output from the walk both, DQW and DCA, is designed to have zero probabilities for an alternate position, however, due to addressing crosstalk in the system, we see a small amount of population in these states. The same mechanism can populate the state \(\left|1000\right\rangle\) of the logical encoding not included in our mapping. In fact, the average experimental population registered in this state is  <2% for the deepest circuits and hence does not affect the results significantly.

Comparison between Dirac kinematics and DCA

We use the explicit time-dependent solution of the one-dimensional Dirac equation provided by Strauch18:

$$\Psi (x,t)=\frac{m{\mathcal{N}}}{\pi }\left(\begin{array}{c}{s}^{-1}{K}_{1}(ms)\left[a+i(t+x)\right]+{K}_{0}(ms)\\ {s}^{-1}{K}_{1}(ms)\left[a+i(t-x)\right]+{K}_{0}(ms)\end{array}\right),$$
(9)

where \(s={[{x}^{2}{(a+it)}^{2}]}^{1/2}\), \({\mathcal{N}}=\sqrt{(\pi /2m)}{[{K}_{1}(2ma)+{K}_{0}(2ma)]}^{-1/2}\) the normalized factor and Kn is the modified Bessel Function of order n, to show the corresponding probability density at time t to the DCA after the time-step t, Fig. 5. The relationship between the mass in the Dirac equation and the coin bias parameter is given by,

$$m\approx \frac{{\theta }_{2}}{1-\frac{{\theta }_{2}^{2}}{2}}.$$
(10)
Fig. 5: Dirac kinematics and DCA.
figure 5

Numerical simulation of the explicit time-dependent solution of the one-dimensional Dirac equation (solid blue) and DCA (yellow bars) at (a) t = 3 and (b) t = 5 with a = 0.4 and θ2 = π/20.

Gate block

The compiler breaks down the gate blocks shown in Fig. 3 (Toffoli-CNOT and Toffoli - Toffoli 4 - CNOT) into native R and XX gates as given by the following circuits, which are optimal in the XX-gate count, Fig. 6. Sketch of the XX-gate is meant to symbolize the two-qubit entangling gate between the outer ions inside a square.

Fig. 6: Gate block.
figure 6

a Toffoli 4--CNOT and b Toffoli 4--CNOT.